% BCA        Use results from MLE for the period 1959:1--2003:4 to
%            (1) plot output, labor and investment for the model and
%            the data and to (2) compute variance decompositions
%            for the 1982 recession. 
%

%            Ellen McGrattan, 2-16-02

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

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

x0=[ ...
  -0.12218441160373
   0.23456126785859
   0.21813418432846
  -1.71044434741722
   0.76890460108055
  -0.02328852467396
  -0.10128969659480
  -0.03063606505969
   0.04705520990984
   0.99352262916894
   0.02531186131043
   0.04233819795798
   0.43178571903678
   0.04171548961226
   1.18316739565328
  -0.00160108354708
  -0.04188834170829
  -0.00452014405740
  -0.01798869654933
   0.99730018412890
   0.01081927622854
  -0.00256942082981
   0.00370721649763
  -0.00500766071435
  -0.00623233345783
  -0.00088560323078
  -0.00505372754194
   0.00194528573172
   0.01479167606597
  -0.00000249657496];

[x1,f,g,code,status]   = uncmin(x0,'mleq',0);
[L,Sbar,P0,P,Q,A,B,C,param] = mleq(x1);
Y0      = 81;
gn      = param(1);
gz      = param(2);
beta    = param(3);
delta   = param(4);
psi     = param(5);
sigma   = param(6);
theta   = param(7);
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);
  Ktp(i,1)   = ((1-delta)*Kt(i)+exp(lxt(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';

lktn     = zeros(T,1);
lktz     = zeros(T,1);
lktl     = zeros(T,1);
lktx     = zeros(T,1);
lktg     = zeros(T,1);
lktn(Y0) = lkt(Y0);
lktz(Y0) = lkt(Y0);
lktl(Y0) = lkt(Y0);
lktx(Y0) = lkt(Y0);
lktg(Y0) = lkt(Y0);
for i=Y0+1:T;
  lktn(i) = A(1,1)*lktn(i-1)+A(1,2:5)*[lzt(Y0),tault(Y0),tauxt(Y0),lgt(Y0)]'+A(1,6);
  lktz(i) = A(1,1)*lktz(i-1)+A(1,2:5)*[lzt(i-1),tault(Y0),tauxt(Y0),lgt(Y0)]'+A(1,6);
  lktl(i) = A(1,1)*lktl(i-1)+A(1,2:5)*[lzt(Y0),tault(i-1),tauxt(Y0),lgt(Y0)]'+A(1,6);
  lktx(i) = A(1,1)*lktx(i-1)+A(1,2:5)*[lzt(Y0),tault(Y0),tauxt(i-1),lgt(Y0)]'+A(1,6);
  lktg(i) = A(1,1)*lktg(i-1)+A(1,2:5)*[lzt(Y0),tault(Y0),tauxt(Y0),lgt(i-1)]'+A(1,6);
end;
for i=Y0-1:-1:1;
  lktn(i) = (lktn(i+1)-A(1,2:5)*[lzt(Y0),tault(Y0),tauxt(Y0),lgt(Y0)]'-A(1,6))/A(1,1);
  lktz(i) = (lktz(i+1)-A(1,2:5)*[lzt(i),tault(Y0),tauxt(Y0),lgt(Y0)]'-A(1,6))/A(1,1);
  lktl(i) = (lktl(i+1)-A(1,2:5)*[lzt(Y0),tault(i),tauxt(Y0),lgt(Y0)]'-A(1,6))/A(1,1);
  lktx(i) = (lktx(i+1)-A(1,2:5)*[lzt(Y0),tault(Y0),tauxt(i),lgt(Y0)]'-A(1,6))/A(1,1);
  lktg(i) = (lktg(i+1)-A(1,2:5)*[lzt(Y0),tault(Y0),tauxt(Y0),lgt(i)]'-A(1,6))/A(1,1);
end;
YMn     = [lktn,ones(T,1)*[lzt(Y0),tault(Y0),tauxt(Y0),lgt(Y0),1]]*C';
YMz     = [lktz-lktn+lkt(Y0),lzt,ones(T,1)*[tault(Y0),tauxt(Y0),lgt(Y0),1]]*C';
YMl     = [lktl-lktn+lkt(Y0),ones(T,1)*lzt(Y0),tault,ones(T,1)*[tauxt(Y0),lgt(Y0),1]]*C';
YMx     = [lktx-lktn+lkt(Y0),ones(T,1)*[lzt(Y0),tault(Y0)],tauxt,ones(T,1)*[lgt(Y0),1]]*C';
YMg     = [lktg-lktn+lkt(Y0),ones(T,1)*[lzt(Y0),tault(Y0),tauxt(Y0)],lgt,ones(T,1)]*C';
YMnox   = YMn+YMz+YMl+YMg-ones(T,1)*[lyt(Y0),lxt(Y0),llt(Y0),lgt(Y0)]*3;
YMnoz   = YMn+YMl+YMx+YMg-ones(T,1)*[lyt(Y0),lxt(Y0),llt(Y0),lgt(Y0)]*3;

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([1973,1985,.9,1.1])

figure(2)
fig2d = [t,[exp(Y(:,1)),exp(Y(:,3)-Y(Y0,3)),exp(Y(:,2))]*100];
fig2mz= [t,[exp(YMz(:,1)),exp(YMz(:,3)-YMz(Y0,3)),exp(YMz(:,2))]*100];
fig2ml= [t,[exp(YMl(:,1)),exp(YMl(:,3)-YMl(Y0,3)),exp(YMl(:,2))]*100];
orient tall
subplot(311)
plot(t,exp(Y(:,1)),t,exp(YMz(:,1)),t,exp(YMl(:,1)))
title('Figure 2. U.S.~Data and Models with an Efficiency or Labor Wedge')
legend('Data','Efficiency wedge','Labor wedge')
axis([1973,1985,.9,1.1])
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([1973,1985,.9,1.1])
subplot(313)
plot(t,exp(Y(:,2)),t,exp(YMz(:,2)),t,exp(YMl(:,2)))
axis([1973,1985,0,.3])

figure(3)
fig3d = [t,[exp(Y(:,1)),exp(Y(:,3)-Y(Y0,3)),exp(Y(:,2))]*100];
fig3mg= [t,[exp(YMg(:,1)),exp(YMg(:,3)-YMg(Y0,3)),exp(YMg(:,2))]*100];
fig3mx= [t,[exp(YMx(:,1)),exp(YMx(:,3)-YMx(Y0,3)),exp(YMx(:,2))]*100];
orient tall
subplot(311)
plot(t,exp(Y(:,1)),t,exp(YMg(:,1)),t,exp(YMx(:,1)))
title('Figure 3. U.S.~Data and Models with an Efficiency or Labor Wedge')
legend('Data','Government wedge','Investment wedge')
axis([1973,1985,.9,1.1])
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([1973,1985,.9,1.1])
subplot(313)
plot(t,exp(Y(:,2)),t,exp(YMg(:,2)),t,exp(YMx(:,2)))
axis([1973,1985,.2,.4])

figure(4)
fig4d   = [t,[exp(Y(:,1)),exp(Y(:,3)-Y(Y0,3)),exp(Y(:,2))]*100];
fig4mnox= [t,[exp(YMnox(:,1)),exp(YMnox(:,3)-YMnox(Y0,3)), ...
              exp(YMnox(:,2))]*100];
fig4mnoz= [t,[exp(YMnoz(:,1)),exp(YMnoz(:,3)-YMnoz(Y0,3)), ...
              exp(YMnoz(:,2))]*100];
orient tall
subplot(311)
plot(t,exp(Y(:,1)),t,exp(YMnox(:,1)),t,exp(YMnoz(:,1)))
title('Figure 4. U.S.~Data and Models without an Investment or Efficiency Wedge')
legend('Data','No investment wedge','No efficiency wedge')
axis([1973,1985,.9,1.1])
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([1973,1985,.9,1.1])
subplot(313)
plot(t,exp(Y(:,2)),t,exp(YMnox(:,2)),t,exp(YMnoz(:,2)))
axis([1973,1985,.2,.4])


%
% Compute variance decompositions and spectra
%
orders     = [2,3,4,5;
              2,3,5,4;
              2,4,3,5;
              2,4,5,3;
              2,5,3,4;
              2,5,4,3;
              3,2,4,5;
              3,2,5,4;
              3,4,2,5;
              3,4,5,2;
              3,5,2,4;
              3,5,4,2;
              4,2,3,5;
              4,2,5,3;
              4,3,2,5;
              4,3,5,2;
              4,5,2,3;
              4,5,3,2;
              5,2,3,4;
              5,2,4,3;
              5,3,2,4;
              5,3,4,2;
              5,4,2,3;
              5,4,3,2];
for norder=1:24;
order      = orders(norder,:);
jz         = find(order==2);
jtaul      = find(order==3);
jtaux      = find(order==4);
jg         = find(order==5);
nsteps     = 50;
nsteps     = 100;
nfreq      = 200;
As         = A([1,order],[1,order]);
Bs         = B([1,order],order-1);
Cs         = C(:,[1,order]);

VX         = dlyap(As,Bs*Bs');
VY         = Cs*VX*Cs';
L          = zeros(5);
L(2:5,2:5) = chol(Bs(2:5,:)*Bs(2:5,:)')';
e2         = zeros(5); e2(2,2)=1;
e3         = zeros(5); e3(3,3)=1;
e4         = zeros(5); e4(4,4)=1;
e5         = zeros(5); e5(5,5)=1;
VX2        = dlyap(As,L*e2*L');
VX3        = dlyap(As,L*e3*L');
VX4        = dlyap(As,L*e4*L');
VX5        = dlyap(As,L*e5*L');
VY2        = Cs*VX2*Cs';
VY3        = Cs*VX3*Cs';
VY4        = Cs*VX4*Cs';
VY5        = Cs*VX5*Cs';

VX         = VX+1e-20;
VY         = VY+1e-20;
VYj        = 0*VY;
VYj2       = 0*VY;
VYj3       = 0*VY;
VYj4       = 0*VY;
VYj5       = 0*VY;
vy_y       = zeros(nsteps,4);
vy_x       = zeros(nsteps,4);
vy_l       = zeros(nsteps,4);
for j=1:nsteps;

  VYj  = VYj  + Cs*As^(j-1)*Bs*Bs'*(As^(j-1))'*Cs';
  VYj2 = VYj2 + Cs*As^(j-1)*L*e2*L'*(As^(j-1))'*Cs';
  VYj3 = VYj3 + Cs*As^(j-1)*L*e3*L'*(As^(j-1))'*Cs';
  VYj4 = VYj4 + Cs*As^(j-1)*L*e4*L'*(As^(j-1))'*Cs';
  VYj5 = VYj5 + Cs*As^(j-1)*L*e5*L'*(As^(j-1))'*Cs';
  tem  = [diag(VYj2./VYj),diag(VYj3./VYj),diag(VYj4./VYj),diag(VYj5./VYj)];
  vy_y(j,:) = tem(1,:);
  vy_x(j,:) = tem(2,:);
  vy_l(j,:) = tem(3,:);

end;

ome = linspace(0.00001,pi,nfreq);
i   = sqrt(-1);
sy_y      = zeros(nfreq,4);
sy_x      = zeros(nfreq,4);
sy_l      = zeros(nfreq,4);
for j=1:length(ome);
  SYj     = Cs*inv(exp(i*ome(j))*eye(5)-As)*Bs*Bs'* ...
               inv(exp(-i*ome(j))*eye(5)-As')*Cs';
  SYj2    = Cs*inv(exp(i*ome(j))*eye(5)-As)*L*e2*L'* ...
               inv(exp(-i*ome(j))*eye(5)-As')*Cs';
  SYj3    = Cs*inv(exp(i*ome(j))*eye(5)-As)*L*e3*L'* ...
               inv(exp(-i*ome(j))*eye(5)-As')*Cs';
  SYj4    = Cs*inv(exp(i*ome(j))*eye(5)-As)*L*e4*L'* ...
               inv(exp(-i*ome(j))*eye(5)-As')*Cs';
  SYj5    = Cs*inv(exp(i*ome(j))*eye(5)-As)*L*e5*L'* ...
               inv(exp(-i*ome(j))*eye(5)-As')*Cs';
  tem     = real([diag(SYj2./SYj),diag(SYj3./SYj),diag(SYj4./SYj), ...
                  diag(SYj5./SYj)]);
  sy_y(j,:) = tem(1,:);
  sy_x(j,:) = tem(2,:);
  sy_l(j,:) = tem(3,:);
end;

syz(:,norder)    = sy_y(:,jz);
sytaul(:,norder) = sy_y(:,jtaul);
sytaux(:,norder) = sy_y(:,jtaux);
syg(:,norder)    = sy_y(:,jg);

end;

figure(5)
fig5 = [[1:nfreq]',ome',[mean(syz')',mean(sytaul')',mean(sytaux')']*100];
orient tall
plot(ome,mean(syz'),ome,mean(sytaul'),ome,mean(sytaux'),ome,mean(syg'))
title('Average Contributions of 4 Shocks to Output Spectrum')
legend('z shock','\tau_l shock','\tau_x shock','g shock')

