% TRUTH2  Compute the true model impulse response function for the 
%         two-shock, unit-root version of the RBC model.  This version
%         assumes that taux and g are constant.
%         The code also computes and plots the traces of B_jB_j' where
%         B_j is the jth VAR coefficient implied by the model.
%       

%         Ellen McGrattan, 3-15-04

%______________________________________________________________________________
% PARAMETERS
%
gn         = (1.01)^(1/4)-1;
gz         = (1.02)^(1/4)-1;
beta       = .98^(1/4);
delta      = 1-(1-.06)^(1/4);
psi        = 2.5;
sigma      = 1.000001;
theta      = .33;
param      = [gn;gz;beta;delta;psi;sigma;theta];
alpha      = 0;
pq         = 1;
if pq==1; pu;  end;  % load Sbar,P,Q for dataset PU
if pq==2; fr1; end;  % load Sbar,P,Q for dataset FR1
if pq==3; fr2; end;  % load Sbar,P,Q for dataset F-R2
if pq==4; cev; end;  % load Sbar,P,Q for dataset CEV
if pq==5; gr;  end;  % load Sbar,P,Q for dataset GR

%______________________________________________________________________________
% COMPUTE EQUILIBRIUM
% a. Steady state
%
tem        = (eye(4)-P)\P0;
zs         = exp(tem(1));
tauls      = tem(2);
tauxs      = tem(3);
gs         = exp(tem(4));
kls        = ((1+tauxs)*(1-beta*zs^(-sigma)*(1-delta))/ ...
              (beta*zs^(-sigma)*theta))^(1/(theta-1))*zs;
A          = kls^(theta-1)*zs^(-theta)-(1+gn)+(1-delta)/zs;
B          = (1-tauls)*(1-theta)*kls^theta*zs^(-theta)/psi;
ks         = (B+gs)/(A+B/kls);
cs         = A*ks-gs;
ls         = ks/kls;
ys         = (ks/zs)^theta*ls^(1-theta);
xs         = ys-cs-gs;
X0s        = [log(ks);log(zs);tauls;tauxs;log(gs);1];
X0         = [X0s;X0s];
Y0         = [log(zs);log(zs);log(ls);log(zs)];

%
% b. Call subroutine with residuals:
%
                                                                                
Z          = [log(ks);log(ks);log(ks);log(zs);log(zs);tauls;tauls;
              tauxs;tauxs;log(gs);log(gs)];
del        = max(abs(Z)*1e-5,1e-8);
for i=1:11;
  Zp       = Z;
  Zm       = Z;
  Zp(i)    = Z(i)+del(i);
  Zm(i)    = Z(i)-del(i);
  dR(i,1)  = (res_wedge_ur(Zp,param)-res_wedge_ur(Zm,param))/(2*del(i));
end;
                                                                                
%
% c. Solution:  log k[t+1] = gamma0 + gammak* log k[t] + gamma* S[t]
%
                                                                                
a0         = dR(1);
a1         = dR(2);
a2         = dR(3);
b0         = dR(4:2:11)';
b1         = dR(5:2:11)';
tem        = roots([a0,a1,a2]);
gammak     = tem(find(abs(tem)<1));
gamma      = -((a0*gammak+a1)*eye(4)+a0*P')\(b0*P+b1)';
gamma0     = (1-gammak)*log(ks)-gamma'*[log(zs);tauls;tauxs;log(gs)];
Gamma      = [gammak;gamma;gamma0];

%
% FINAL SOLUTION:  log k[t+1] = A0 + A* log k[t] + B* S[t]
%
%______________________________________________________________________________
                                                                                
%______________________________________________________________________________
%
% EQUATIONS FOR VAR:

gamz       = gamma(1);
gaml       = gamma(2);
philh      =-(psi*ys*(1-theta)+(1-theta)*(1-tauls)*ys*(1-ls)/ls*theta+ ...
             (1-theta)*(1-tauls)*ys);
philk      = (psi*ys*theta+psi*(1-delta)*ks/zs- ...
             (1-theta)*(1-tauls)*ys*(1-ls)/ls *theta)/philh;
philz      = (psi*ys*(-theta)-psi*(1-delta)*ks/zs+ ...
             (1-theta)*theta*(1-tauls)*ys*(1-ls)/ls)/philh;
phill      = ((1-theta)*(1-tauls)*ys*(1-ls)/ls *(1/(1-tauls)))/philh;
philg      = (-psi*gs)/philh;
philkp     = (-psi*(1+gn)*ks)/philh;
psiz       = philz+philkp*gamz;
psil       = phill+philkp*gaml;

D0         = [1-theta-theta*(philz+philkp*gamz), ...
                          -theta*(phill+philkp*gaml)
              philz+philkp*gamz,phill+philkp*gaml];
D1         = [theta*(1+philz-(philk-1)*gamz-philkp*gamz*(gammak-1)), ...
                   theta*(phill-(philk-1)*gaml-philkp*gaml*(gammak-1));
              -(alpha*philz-philk*gamz-philkp*gamz*(gammak-alpha)), ...
                  -(alpha*phill-philk*gaml-philkp*gaml*(gammak-alpha))];
D2         = [-theta*gamz*(gammak-1)*(philkp*gammak+philk-1), ...
                -theta*gaml*(gammak-1)*(philkp*gammak+philk-1);
             gamz*(gammak-alpha)*(philkp*gammak+philk), ...
                 gaml*(gammak-alpha)*(philkp*gammak+philk)];
D3         = D2*gammak;
D4         = D3*gammak;
D5         = D4*gammak;
P          = P(1:2,1:2);
Q          = Q(1:2,1:2);
C1         = (D0*P+D1)*inv(D0);
C2         = (D0*P^2+D1*P+D2)*inv(D0);
C3         = C2*D0*P*inv(D0)+D3*inv(D0);
B1         = C1;
B2         = C2-B1*C1;
B3         = C3-B1*C2-B2*C1;
Tr         = [sqrt(trace(B1*B1'));
              sqrt(trace(B2*B2'));
              sqrt(trace(B3*B3'))];
S2         = cumsum([sqrt(2);
              sqrt(trace(B1*B1'));
              sqrt(trace(B2*B2'));
              sqrt(trace(B3*B3'))]);
sum2       = S2(4);
S2         = S2(2:4);
ntraces    = 100;
Cbar0      = eye(2)+C1+C2+C3;
for tt=4:ntraces;
  tts = int2str(tt);
  tts1= int2str(tt-1);
  eval(['D',tts,'=gammak*D',tts1,';'])
  eval(['C',tts,'=C',tts1,'*D0*P*inv(D0)+D',tts,'*inv(D0);'])
  sum1=0;
  for ii=1:tt-1; 
    iis = int2str(ii);
    its = int2str(tt-ii);
    eval(['sum1 = sum1 + B',iis,'*C',its,';'])
  end; 
  eval(['B',tts,'=C',tts,'-sum1;']);
  eval(['sum2 = sum2+sqrt(trace(B',tts,'*B',tts,'''));'])
  eval(['Tr(tt) = sqrt(trace(B',tts,'*B',tts,'''));'])
  S2(tt)=sum2;
  eval(['tem1=C',tts,'*D0*Q;'])
  eval(['Cbar0=Cbar0+C',tts,';'])
end; 

A0         = D0*Q;
A1         = D1*Q;
A2         = D2*Q;
A3         = D3*Q;
AA         = [A0(:,1),A1(:,1),A2(:,1)*[1,gammak.^[1:98]]]';
figure(1)
for ii=1:length(AA);
   AA2(ii,1) = (alpha.^[ii-1:-1:0])*AA(1:ii,2);
end;
IR         = [cumsum(AA(:,1)),AA2]*100;
plot(0:100,IR);
title('Impulse responses to shock in technology')
legend('Productivity','Hours')
figure(2)
plot(S2(2:ntraces)-S2(1:ntraces-1))
title('Change in the sum of the traces')


