% INITPW     Compute MLE estimates for the period 1901--2000 to get
%            initial estimates for the postwar estimation.
%

%            Ellen McGrattan, 2-16-02
%            Revised, ERM, 11-8-06

%------------------------------------------------------------------------------
load uszvar1.dat
global ZVAR
nt   = input('Number of periods (40 or 100)? ');
t    = uszvar1(1:nt,1);
ZVAR = uszvar1(1:nt,2:5);
x0a=[ ...
   0.54082135494385   0.53688220711471   0.53632298081285
  -0.18961520211220  -0.19633927001699  -0.20039531347703
   0.28589191522520   0.29058999479661   0.29210706768606
  -2.79358307142896  -2.79623988122216  -2.80780348317646
   0.73194722119709   0.56390392825651   0.43182517280778
  -0.14950077057243   0.04695540205288   0.17815955068537
  -0.01135154662368  -0.37800322911024  -0.65916167358026
   0.05206442487424   0.08977302569669   0.12063917741983
   1.03867122470732   0.99533700152441   0.96304256251874
  -0.01965637369407   0.05448227333478   0.11170155154008
  -0.31698595387468  -0.19009624055368  -0.08657508959716
   0.39035390644605   0.21919090712206   0.09382053325445
   0.07313270008334   0.33666352061755   0.55959882698031
   0.74982873714821   0.76639431904874   0.78311669515450
  -0.05746385259960  -0.05675342255460   0.05614910291263
   0.00561454650003   0.00425205357843  -0.00326537888522
  -0.00298795420604   0.04319778764689  -0.17603674722045
   0.05551125567149   0.05542040154641   0.05526762622885
  -0.00025361856807   0.02058061530191   0.08107289332405
  -0.03694910344636  -0.07681637702580  -0.18875619010229
   0.22143490701858   0.22112658504395   0.22088105811702]; % t=1901:1940
%  adja = 0                3.22               12.88

x0b=[ ...
   0.56129229542991   0.53361937811488   0.53845530645040
  -0.22315683391720  -0.18731348652782  -0.18358932451124
   0.35940299493795   0.26493094444701   0.27053055037658
  -2.31053250090644  -2.69213183271747  -2.77729399525531
   0.90684739795147   0.92022994446255   0.93481145145186
  -0.03122926159168   0.01761941263284   0.06311171605897
  -0.02108604414666  -0.03910923936791  -0.15357502558886
   0.06759103409288   0.06173502693479   0.05094880731915
   1.06064502909895   1.00624450123758   0.97010706688692
  -0.03635578467867  -0.00801280802084   0.03514397572032
   0.05474673522732   0.04349311619849   0.01821218949912
   0.17843437665379   0.13800062069869   0.06926955105946
   0.82390093717686   0.78918690892135   0.74307481445789
   0.87544873454186   0.94821840408492   0.96444356405522
  -0.04593755989777  -0.04597420974653  -0.04619158158676
   0.01967151596379   0.01917525092030   0.01870284601253
  -0.02478501995215  -0.01236378739766   0.02647798195864
   0.05395834709159   0.05197393897980   0.05143250999846
  -0.03207776032144  -0.05120198182181  -0.08591038653482
  -0.05305289195720  -0.09709413673436  -0.23289864958749
   0.21223072021859   0.19991098002854   0.19800381242284]; % t=1901:2000
%  adja = 0                3.22               12.88


adja = input('Annual adjustment cost parameter (0, 3.22, or 12.88) = ');
if (adja==0);
  icol = 1;
elseif adja==3.22;
  icol = 2;
else 
  icol = 3;
end;
if (nt==40);
  x0 = x0a(:,icol);
else
  x0 = x0b(:,icol);
end;

[x,f,g,code,status]  = uncmin(x0,'mle1adj',adja);
[L,Sbar,P0,P,Q,A,B,C,param] = mle1adj(x,adja);
                                                                                
[P0q,Pq,Qq]=ara2arq(P0,P,Q);
