% PWBCA      Use results from MLE for the period 1959:1--2004:3 to
%            (1) plot output, labor and investment for the model and
%            the data and to (2) compute variance decompositions
%            for the 1982 recession.  This code also allows for 
%            3 levels of adjustment costs
%

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

%------------------------------------------------------------------------------

load uszvarq.dat
global ZVAR
t    = uszvarq(:,1);
ZVAR = uszvarq(:,2:5);

x0=[ ...
  -0.02392130483182  -0.02593003143300  -0.03868373032989
   0.32793968650403   0.32528078240235   0.32533115894097
   0.48344057752536   0.47640296988723   0.42386985342091
  -1.53442333795160  -1.52572392233495  -1.53913101241965
   0.97998737588446   0.98814714427224   1.01420535776153
  -0.03298062237437  -0.00393620952383   0.00700407325041
  -0.07024526786903   0.00739476978480  -0.10552749200191
   0.00481049463613  -0.00154763964840  -0.03794064209269
  -0.01378479653021   0.02347540555172  -0.04443387632167
   0.95638310701915   0.99522631431640   0.95319275816302
  -0.04600503838943  -0.07813326483521   0.10009893037604
  -0.00810596308553  -0.01292726594313  -0.01576973857835
  -0.01172679028012   0.05213951101474   0.03612391751414
  -0.04508411234601   0.03785132235434   0.01732329566843
   0.89618870904417   0.87933464895897   0.82961577117656
   0.04883996944559   0.04800537603439  -0.03018624159039
   0.01923840639928  -0.05623786887083  -0.01684629253098
   0.05690492926449  -0.02974898581502   0.00403375819985
   0.10407546546790   0.14558343713444   0.13978989734296
   0.97107670610686   0.96500911999557   1.04444661366546
   0.01161970401808   0.01143591934237   0.01137088090910
   0.00141164823025   0.00132263476266   0.00148119047057
  -0.01049727181035  -0.00924790891037  -0.02717953211122
  -0.00057540165855  -0.00038806619786  -0.00002984273316
   0.00644004245925   0.00640541868455   0.00626157343799
   0.00103165759873   0.00264329092744   0.01220319595520
   0.00611247597734   0.00619191723283   0.00589079806760
   0.01584164042916   0.02206456706386   0.02431630028343
   0.01417545102724   0.01348046236832   0.00929368330591
   0.00458359733815   0.00697830381845   0.01154144364817];
Adja = [0,4*3.22,4*12.88];

iadj = input('No adjustment costs (1), BGG level (2), or 4 x BGG (3)? ');
adja = Adja(iadj);
x0   = x0(:,iadj);
[L,Sbar,P0,P,Q,A,B,C,param] = mleqadj(x0,adja);
Y0       = 81;
gn       = param(1);
gz       = param(2);
beta     = param(3);
delta    = param(4);
psi      = param(5);
sigma    = param(6);
theta    = param(7);
adja     = param(8);
adjb     = (1+gn)*(1+gz)-1+delta;
z        = exp(Sbar(1));
taul     = Sbar(2);
taux     = Sbar(3);
g        = exp(Sbar(4));
betah    = beta*(1+gz)^(-sigma);
kl       = ((1+taux)*(1-betah*(1-delta))/(betah*theta))^(1/(theta-1))*z;
yk       = (kl/z)^(theta-1);
xi1      = yk-(1+gz)*(1+gn)+1-delta;
xi2      = (1-taul)*(1-theta)*(kl)^theta*z^(1-theta)/psi;
k        = (xi2+g)/(xi1+xi2/kl);
c        = xi1*k-g;
l        = k/kl;
y        = yk*k;
x        = y-c-g;
lk       = log(k);
lc       = log(c);
ll       = log(l);
ly       = log(y);
lx       = log(x);
lg       = log(g);
lz       = log(z);
T        = length(ZVAR);
Y        = log(ZVAR)-log([(1+gz).^[0:T-1]',(1+gz).^[0:T-1]', ...
                         ones(T,1),(1+gz).^[0:T-1]']);
lyt      = Y(:,1);
lxt      = Y(:,2);
llt      = Y(:,3);
lgt      = Y(:,4);
lkt(1,1) = lk;
Kt(1,1)  = exp(lk);

for i=1:T;
  lktp(i,1)  = lk+((1-delta)*(lkt(i)-lk)+x/k*(lxt(i)-lx))/(1+gz)/(1+gn);
  lkt(i+1,1) = lktp(i);
  phit(i,1)  = adja/2*(exp(lxt(i))/Kt(i)-adjb)^2;
  Ktp(i,1)   = ((1-delta)*Kt(i)+exp(lxt(i))-phit(i)*Kt(i))/(1+gz)/(1+gn);
  Kt(i+1,1)  = Ktp(i);
end;
lkt      = lkt(1:T);
Kt       = Kt(1:T);
lct      = lc+(y*(lyt-ly)-x*(lxt-lx)-g*(lgt-lg))/c;
lzt      = lz+(lyt-ly-theta*(lkt-lk))/(1-theta)-llt+ll;
tault    = taul+(1-taul)*(lyt-ly-lct+lc-1/(1-l)*(llt-ll));
tauxt    = (lxt-C(2,1)*lkt-C(2,2)*lzt-C(2,3)*tault-C(2,5)*lgt-C(2,6))/C(2,4);
tauxchk  = (lyt-C(1,1)*lkt-C(1,2)*lzt-C(1,3)*tault-C(1,5)*lgt-C(1,6))/C(1,4);
Ct       = exp(lyt)-exp(lxt)-exp(lgt);
Zt       = (exp(lyt)./(Kt.^theta.*exp(llt).^(1-theta))).^(1/(1-theta));
Tault    = 1-psi/(1-theta)* (Ct./exp(lyt)) .*(exp(llt)./(1-exp(llt)));
Xt0      = [lkt,lzt,tault,tauxt,lgt,ones(T,1)];
YM0      = Xt0*C';

s0       = [lzt(Y0);tault(Y0);tauxt(Y0);lgt(Y0)];

[i,j,C0] = fixexpadj(Sbar,P,Q,s0,[0;0;0;0],adja);
[i,j,C1] = fixexpadj(Sbar,P,Q,s0,[1;0;0;0],adja);
[i,j,C2] = fixexpadj(Sbar,P,Q,s0,[0;1;0;0],adja);
[i,j,C3] = fixexpadj(Sbar,P,Q,s0,[0;0;1;0],adja);
[i,j,C4] = fixexpadj(Sbar,P,Q,s0,[0;0;0;1],adja);
o        = ones(183,1);
YMn      = (Xt0-o*Xt0(Y0,:))*C0'+o*YM0(Y0,:);
YMz      = (Xt0-o*Xt0(Y0,:))*(C1-C0)'+o*YM0(Y0,:);
YMl      = (Xt0-o*Xt0(Y0,:))*(C2-C0)'+o*YM0(Y0,:);
YMx      = (Xt0-o*Xt0(Y0,:))*(C3-C0)'+o*YM0(Y0,:);
YMg      = (Xt0-o*Xt0(Y0,:))*(C4-C0)'+o*YM0(Y0,:);
YMall    = (Xt0-o*Xt0(Y0,:))*C'+o*YM0(Y0,:);
YMnox    = (Xt0-o*Xt0(Y0,:))*(C1+C2+C4-2*C0)'+o*YM0(Y0,:);
YMnoz    = (Xt0-o*Xt0(Y0,:))*(C2+C3+C4-2*C0)'+o*YM0(Y0,:);

figure(1)
fig1 =[t,[exp(lyt-lyt(Y0)),(Zt/Zt(Y0)).^(1-theta), ...
       (1-Tault)/(1-Tault(Y0)), ...
       (1+tauxt(Y0))*ones(T,1)./(1+tauxt)]*100];
orient portrait
plot(t,exp(lyt-lyt(Y0)), ...
     t,(Zt/Zt(Y0)).^(1-theta), ...
     t,(1-Tault)/(1-Tault(Y0)), ...
     t,(1+tauxt(Y0))*ones(T,1)./(1+tauxt))
title('Figure 1. U.S. Output and Measured Wedges')
legend('Data','Efficiency Wedge','Labor Wedge','Investment Wedge')
axis([1979,1986,.8,1.2])

figure(2)
fig2d = [t,[exp(Y(:,1)-Y(Y0,1)),exp(Y(:,3)-Y(Y0,3)),exp(Y(:,2))]*100];
fig2mz= [t,[exp(YMz(:,1)-YMz(Y0,1)),exp(YMz(:,3)-YMz(Y0,3)),exp(YMz(:,2))]*100];
fig2ml= [t,[exp(YMl(:,1)-YMl(Y0,1)),exp(YMl(:,3)-YMl(Y0,3)),exp(YMl(:,2))]*100];
orient tall
subplot(311)
plot(t,exp(Y(:,1)-Y(Y0,1)),t,exp(YMz(:,1)-YMz(Y0,1)),t,exp(YMl(:,1)-YMl(Y0,1)))
title('Figure 2. U.S.~Data and Models with an Efficiency or Labor Wedge')
legend('Data','Efficiency wedge','Labor wedge')
axis([1979,1986,.8,1.2])
subplot(312)
plot(t,exp(Y(:,3)-Y(Y0,3)),t,exp(YMz(:,3)-YMz(Y0,3)),t,exp(YMl(:,3)-YMl(Y0,3)))
axis([1979,1986,.9,1.1])
subplot(313)
plot(t,exp(Y(:,2)),t,exp(YMz(:,2)),t,exp(YMl(:,2)))
axis([1979,1986,.1,.4])

figure(3)
fig3d = [t,[exp(Y(:,1)-Y(Y0,1)),exp(Y(:,3)-Y(Y0,3)),exp(Y(:,2))]*100];
fig3mg= [t,[exp(YMg(:,1)-YMg(Y0,1)),exp(YMg(:,3)-YMg(Y0,3)),exp(YMg(:,2))]*100];
fig3mx= [t,[exp(YMx(:,1)-YMx(Y0,1)),exp(YMx(:,3)-YMx(Y0,3)),exp(YMx(:,2))]*100];
orient tall
subplot(311)
plot(t,exp(Y(:,1)-Y(Y0,1)),t,exp(YMg(:,1)-YMg(Y0,1)),t,exp(YMx(:,1)-YMx(Y0,1)))
title('Figure 3. U.S.~Data and Models with an Investment or Government Wedge')
legend('Data','Government wedge','Investment wedge')
axis([1979,1986,.8,1.2])
subplot(312)
plot(t,exp(Y(:,3)-Y(Y0,3)),t,exp(YMg(:,3)-YMg(Y0,3)),t,exp(YMx(:,3)-YMx(Y0,3)))
axis([1979,1986,.8,1.2])
subplot(313)
plot(t,exp(Y(:,2)),t,exp(YMg(:,2)),t,exp(YMx(:,2)))
axis([1979,1986,.1,.4])

figure(4)
fig4d   = [t,[exp(Y(:,1)-Y(Y0,1)),exp(Y(:,3)-Y(Y0,3)),exp(Y(:,2))]*100];
fig4mnox= [t,[exp(YMnox(:,1)-YMnox(Y0,1)),exp(YMnox(:,3)-YMnox(Y0,3)), ...
              exp(YMnox(:,2))]*100];
fig4mnoz= [t,[exp(YMnoz(:,1)-YMnoz(Y0,1)),exp(YMnoz(:,3)-YMnoz(Y0,3)), ...
              exp(YMnoz(:,2))]*100];
orient tall
subplot(311)
plot(t,exp(Y(:,1)-Y(Y0,1)),t,exp(YMnox(:,1)-YMnox(Y0,1)), ...
     t,exp(YMnoz(:,1)-YMnoz(Y0,1)))
title(['Figure 4. U.S.~Data and Models without an Investment or ', ...
                  'Efficiency Wedge'])
legend('Data','No investment wedge','No efficiency wedge')
axis([1979,1986,.8,1.2])
subplot(312)
plot(t,exp(Y(:,3)-Y(Y0,3)),t,exp(YMnox(:,3)-YMnox(Y0,3)), ...
     t,exp(YMnoz(:,3)-YMnoz(Y0,3)))
axis([1979,1986,.8,1.2])
subplot(313)
plot(t,exp(Y(:,2)),t,exp(YMnox(:,2)),t,exp(YMnoz(:,2)))
axis([1979,1986,.1,.4])


