%GMODEL2    computes equilibria for a growth model with two sectors --
%           a corporate sector and a noncorporate sector.  The difference
%           between this version and gmodel.m is the specification of 
%           technology.  Here it is an AR(1).
%
%           Households solve:
%
%               max  E [sum  beta^t U(c[t],l[t])
%              c,l,i
%                   s.t. c[t]+i_c[t]+i_n[t] 
%                            = (1-tau_c[t]) r_c[t]k_c[t]+tau_c[t] delc k_c[t]
%                             +(1-tau_n[t]) r_n[t]k_n[t]+tau_n[t] deln k_n[t]
%                             +(1-tau_h[t]) w[t] h[t] + T[t]
%                        k_c[t+1] = (1-delc) k_c[t] + i_c[t]
%                        k_n[t+1] = (1-deln) k_n[t] + i_n[t]
%                        h[t]+l[t]= 1 
%
%           Corporate firms solve:
%
%               max   p_c[t] y_c[t] -r_c[t] k_c[t] - w[t] h_c[t]   
%             k_c,h_c  
%                     s.t. y[t]=k_c[t]^theta[t] (z[t] mu^t h_c[t])^(1-theta[t])
%
%
%           Noncorporate firms solve:
%
%               max   p_n[t] y_n[t] -r_n[t] k_n[t] - w[t] h_n[t]
%             k_n,h_n  
%                     s.t. y[t]=k_n[t]^phi[t] (z[t] mu^t h_n[t])^(1-phi[t])
%
%           Additionally,
%
%               y[t] = c[t] + i_c[t] + i_n[t] + g[t]
%               h[t] = h_c[t] + h_n[t]
%               g[t] + T[t] = tau_c[t] (r_c[t]-delc) k_c[t]
%                            +tau_n[t] (r_n[t]-deln) k_n[t] + tau_h[t]w[t]h[t]
%
%           Stochastic variables:  
%               ln z[t+1]  = rho_z ln z[t] + eta[t+1]
%               tau_c[t+1] = (1-rho_c) mtauc  + rho_c tau_c[t] + eps_c[t+1]
%               tau_n[t+1] = (1-rho_n) mtaun  + rho_n tau_n[t] + eps_n[t+1]
%               tau_h[t+1] = (1-rho_h) mtauh  + rho_h tau_h[t] + eps_h[t+1]
%               ln g[t+1]  = (1-rho_g) mg     + rho_g ln g[t]  + eps_g[t+1]
%               theta[t+1] = (1-rho_t) mtheta + rho_t theta[t] + eps_t[t+1]
%               phi[t+1]   = (1-rho_p) mphi   + rho_p phi[t]   + eps_p[t+1]
%
%           Computation: find h_c(k_c,k_n,s), i_c(k_c,k_n,s)
%                        where s = [z,tau_c,tau_n,tau_h,ln g,theta,phi]'
%                        NOTE: all nonstationary variables must be divided
%                              by mu^t before computing. 
%             
%

%           Ellen McGrattan, 5-22-00
%           Revised, ERM, 5-22-00

%______________________________________________________________________________
%
%PARAMETERS:
%
%  Preferences:
%
%    U(c,l)   = (c l^psi)^(1-sigma)/(1-sigma)
%    beta     = the pre-growth discount factor 
%   
%  Technology:
%
%    y_c      = k_c^theta * (z h_c)^(1-theta)       = corporate production
%    y_n      = k_n^phi   * (z h_n)^(1-phi)         = noncorporate production
%    y        = (a1 y_c^rho + a2 y_n^rho)^(1/rho)   = GNP
%
%  Capital accumulation
%
%    k_c[t+1] = [(1-delc) k_c[t] + i_c[t]]/mu
%    k_n[t+1] = [(1-deln) k_n[t] + i_n[t]]/mu
%
%______________________________________________________________________________
%

clear Re

psi      = .01;
sigma    = 5;
beta     = .99*(1.0045)^(sigma-1);
a1       = .5;
a2       = .5;
rho      = .5;
delc     = .025;
deln     = .025;
mu       = 1.005;

mtauc    = .0;
mtaun    = .0; 
mtauh    = .0;
mg       = .0001;
mthet    = .36;
mphi     = .4;

rhoz     = .9992;
rhoz     = .95;
rhoc     = .0;
rhon     = .0;
rhoh     = .0;
rhog     = .0;
rhot     = .0;
rhop     = .0;

sigz     = .02;
sigc     = 0.0;
sign     = 0.0;
sigh     = 0.0;
sigg     = 0.0;
sigt     = 0.0;
sigp     = 0.0;
cortp    = 0;

param    = [a1;a2;beta;delc;deln;mphi;mthet;mu;psi;rho;sigma];
betm     = beta*mu^(-sigma);
mlg      = log(mg);

%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

ycy    = .6;   % Guess of corporate GDP to total GDP
hc     = .4;  % Guess of hours of work in the corporate sector

z0     = [ycy;hc];
del    = diag(z0) * .00001;

for i=1:20;

   ycy = z0(1); 
   hc  = z0(2); 
   yny = ((1-a1*ycy^rho)/a2)^(1/rho);
   pc  = a1*ycy^(rho-1);
   pn  = a2*yny^(rho-1);
   hkc = ((1-betm+(1-mtauc)*betm*delc)/(betm*(1-mtauc)*mthet*pc))^(1/(1-mthet));
   hkn = ((1-betm+(1-mtaun)*betm*deln)/(betm*(1-mtaun)*mphi*pn))^(1/(1-mphi));
   kc  = hc/hkc;
   yc  = hkc^(1-mthet)*kc;
   y   = yc/ycy;
   yn  = (y^rho/a2-a1/a2*yc^rho)^(1/rho);
   kn  = yn*hkn^(mphi-1);
   hn  = hkn*kn; 
   ic  = (mu-1+delc)*kc;
   in  = (mu-1+deln)*kn;
   c   = y-ic-in-mg;
    
   f   = [pc*(1-mthet)*hkc^(-mthet)-pn*(1-mphi)*hkn^(-mphi);
          psi*c-(1-mtauh)*pc*(1-mthet)*hkc^(-mthet)*(1-hc-hn)];

   ycy = z0(1) + del(1,1);
   yny = ((1-a1*ycy^rho)/a2)^(1/rho);
   pc  = a1*ycy^(rho-1);
   pn  = a2*yny^(rho-1);
   hkc = ((1-betm+(1-mtauc)*betm*delc)/(betm*(1-mtauc)*mthet*pc))^(1/(1-mthet));
   hkn = ((1-betm+(1-mtaun)*betm*deln)/(betm*(1-mtaun)*mphi*pn))^(1/(1-mphi));
   kc  = hc/hkc;
   yc  = hkc^(1-mthet)*kc;
   y   = yc/ycy;
   yn  = (y^rho/a2-a1/a2*yc^rho)^(1/rho);
   kn  = yn*hkn^(mphi-1);
   hn  = hkn*kn; 
   ic  = (mu-1+delc)*kc;
   in  = (mu-1+deln)*kn;
   c   = y-ic-in-mg;
    
   df(:,1) ...
       = ([pc*(1-mthet)*hkc^(-mthet)-pn*(1-mphi)*hkn^(-mphi);
           psi*c-(1-mtauh)*pc*(1-mthet)*hkc^(-mthet)*(1-hc-hn)]-f)/del(1,1);

   ycy = z0(1);
   hc  = z0(2) + del(2,2);
   yny = ((1-a1*ycy^rho)/a2)^(1/rho);
   pc  = a1*ycy^(rho-1);
   pn  = a2*yny^(rho-1);
   hkc = ((1-betm+(1-mtauc)*betm*delc)/(betm*(1-mtauc)*mthet*pc))^(1/(1-mthet));
   hkn = ((1-betm+(1-mtaun)*betm*deln)/(betm*(1-mtaun)*mphi*pn))^(1/(1-mphi));
   kc  = hc/hkc;
   yc  = hkc^(1-mthet)*kc;
   y   = yc/ycy;
   yn  = (y^rho/a2-a1/a2*yc^rho)^(1/rho);
   kn  = yn*hkn^(mphi-1);
   hn  = hkn*kn; 
   ic  = (mu-1+delc)*kc;
   in  = (mu-1+deln)*kn;
   c   = y-ic-in-mg;
    
   df(:,2) ...
       = ([pc*(1-mthet)*hkc^(-mthet)-pn*(1-mphi)*hkn^(-mphi);
           psi*c-(1-mtauh)*pc*(1-mthet)*hkc^(-mthet)*(1-hc-hn)]-f)/del(2,2);

   z0  = z0 - df\f;

end;
ycy = z0(1);
hc  = z0(2);
yny = ((1-a1*ycy^rho)/a2)^(1/rho);
pc  = a1*ycy^(rho-1);
pn  = a2*yny^(rho-1);
hkc = ((1-betm+(1-mtauc)*betm*delc)/(betm*(1-mtauc)*mthet*pc))^(1/(1-mthet));
hkn = ((1-betm+(1-mtaun)*betm*deln)/(betm*(1-mtaun)*mphi*pn))^(1/(1-mphi));
kc  = hc/hkc;
yc  = hkc^(1-mthet)*kc;
y   = yc/ycy;
yn  = (y^rho/a2-a1/a2*yc^rho)^(1/rho);
kn  = yn*hkn^(mphi-1);
hn  = hkn*kn; 
ic  = (mu-1+delc)*kc;
in  = (mu-1+deln)*kn;
c   = y-ic-in-mg;
rf  = mu^sigma/beta-1;
rc  = pc*mthet*yc/kc;

%______________________________________________________________________________




%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% The ordering of variables to the residuals is
%
%     Z  =  [ k_{c,t+1}, k_{n,t+1}, h_{c,t+1}, i_{c,t+1}, 
%               k_{c,t}, k_{n,t}, h_{c,t}, i_{c,t},
%               eta_{t+1}, tau_{c,t+1}, tau_{n,t+1}, tau_{h,t+1}, g_{t+1},
%               theta_{t+1}, phi_{t+1}, eta_t, tau_{c,t}, tau_{n,t}, 
%               tau_{h,t}, g_t, theta_t, phi_t]'
%
%

Z        = [log([kc,kn,hc,ic,kc,kn,hc,ic]), ...
            0,mtauc,mtaun,mtauh,mlg,mthet,mphi, ...
            0,mtauc,mtaun,mtauh,mlg,mthet,mphi]';
nZ       = 22;
del      = max(abs(Z)*1e-4,1e-8);
disp(sprintf(['Calculating ',int2str(nZ),' numerical derivatives ...'],nZ))
for i=1:nZ;
  Zp       = Z;
  Zm       = Z;
  Zp(i)    = Zp(i)+del(i);
  Zm(i)    = Zm(i)-del(i);
  [Rp,R2p] = res_gmodel2(Zp,param);
  [Rm,R2m] = res_gmodel2(Zm,param);
  dR(:,i)  = (Rp-Rm)/(2*del(i));
  dR2(:,i) = (R2p-R2m)/(2*del(i));
end;
disp('Computing equilibrium ...')
%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF DETERMINISTIC PART OF PROBLEM (A)
% 

%
% The residuals are given by
%
%    E_t { Psi1 X_{t+1} + Psi2 D_{t+1} + Psi3 X_t + Psi4 D_t 
%                                      + Psi5 S_{t+1} + Psi6 S_t } = 0
%
%          Phi1 X_{t+1} + Phi2 D_{t+1} + Phi3 X_t + Phi4 D_t 
%                                      + Phi5 S_{t+1} + Phi6 S_t   = 0
%
% where 
%        X_t = [ k_{c,t} ]   D_t = [ h_{c,t} ]   S_t = [ eta_t     ]
%              [ k_{n,t} ]         [ i_{c,t} ]         [ tau_{c,t} ]
%                                                      [ tau_{n,t} ]
%                                                      [ tau_{h,t} ]
%                                                      [ g_t       ]
%                                                      [ theta_t   ]
%                                                      [ phi_t     ]
%
%
% The solution to the deterministic part of the problem is given by
%
%     D_t = A X_t
%
%

Psi1  = dR(1:2,1:2);
Psi2  = dR(1:2,3:4);
Psi3  = dR(1:2,5:6);
Psi4  = dR(1:2,7:8);
Psi5  = dR(1:2,9:15);
Psi6  = dR(1:2,16:22);

Phi1  = dR(3:4,1:2);
Phi2  = dR(3:4,3:4);
Phi3  = dR(3:4,5:6);
Phi4  = dR(3:4,7:8);
Phi5  = dR(3:4,9:15);
Phi6  = dR(3:4,16:22);

alp1  = dR2(5:6);
alp2  = dR2(7:8);
alp3  = dR2(16:22);

[V,L] = eig(-[Phi3,Phi4;Psi3,Psi4],[Phi1,Phi2;Psi1,Psi2]);
L     = diag(L);
[D,I] = sort(abs(L));
if (~( (D(2) <= 1) & (D(3) > 1) ) )
  i      = sum(D<=1);
  disp(sprintf( ...
  'WARNING: program assumes %g eigenvalues inside unit circle, found %g',2,i));
end;
V11   = V(1:2,I(1:2));
V21   = V(3:4,I(1:2));
A     = real(V21/V11);


%______________________________________________________________________________


%______________________________________________________________________________
%
% SOLUTION OF COEFFICIENTS ON STOCHASTIC VARIABLES (B)
%

P     = [rhoz,0,0,0,0,0,0;
         0,rhoc,0,0,0,0,0;
         0,0,rhon,0,0,0,0;
         0,0,0,rhoh,0,0,0;
         0,0,0,0,rhog,0,0;
         0,0,0,0,0,rhot,0;
         0,0,0,0,0,0,rhop];
vB    = (kron(P',Psi2)+kron(eye(7),-(Psi1+Psi2*A)*(Phi1\Phi4))+ ...
          kron(eye(7),Psi4))\vec((Psi1+Psi2*A)*(Phi1\Phi5)*P-Psi5*P+ ...
                                 (Psi1+Psi2*A)*(Phi1\Phi6)-Psi6);

B     = zeros(2,7);
B(:)  = vB;

%
% FINAL SOLUTION:  
%     D[t]   = A* X[t] + B* S[t]
%     S[t+1] = P* S[t] + Sig_e eps[t+1]  
%     X[t+1] = -Phi1\[ (Phi3+Phi4 A) X[t] +Phi5 S[t+1] +(Phi6+Phi4 B) S[t] ]
%            = -Phi1\[ (Phi3+Phi4 A) X[t] +(Phi5 P+Phi6+Phi4 B) S[t] 
%                                         + Phi5 Sig_e eps[t+1] ]
%______________________________________________________________________________


%______________________________________________________________________________
% 
% IMPULSE RESPONSES FOR SYSTEM
%
%       [ X[t+1] ] = [ Ao  Bo Means]  [ X[t] ]      [Co   ] 
%       [ S[t+1] ]   [ 0   P  Means]  [ S[t] ]   +  [Sig_e] * e[t+1]
%       [   1    ]   [ 0   0   1   ]  [   1  ]      [0    ]
%
%                         (G)                        (H)         

disp('Generating impulse responses ...')

Sige  = [sigz,0,0,0,0,0,0;
         0,sigc,0,0,0,0,0;
         0,0,sign,0,0,0,0;
         0,0,0,sigh,0,0,0;
         0,0,0,0,sigg,0,0;
         0,0,0,0,0,sigt,cortp*sigt/sqrt(1-cortp*cortp);
         0,0,0,0,0,0,sigp/sqrt(1-cortp*cortp)];
G     = [-Phi1\(Phi3+Phi4*A),-Phi1\(Phi5*P+Phi6+Phi4*B);
          zeros(7,2), P];
H     = [ -(Phi1\Phi5)* Sige;
          Sige
          zeros(1,7)];
M     = [log(kc);log(kn);0;mtauc;mtaun;mtauh;mlg;mthet;mphi];
G     = [G, (eye(9)-G)*M;
         zeros(1,9),1];


alp   = [alp1+alp2*A,alp2*B+alp3]';
e     = zeros(9,1);
e(3)  = 1;
gam   = (alp'*(eye(9)-G(1:9,1:9)))';
Vx    = dlyap(G(1:9,1:9),H(1:9,:)*H(1:9,:)');
Erf   = exp((gam'*Vx*gam-alp'*H(1:9,:)*H(1:9,:)'*alp)/2)*mu^sigma/beta-1;


randn('seed',10462)
nsim  = 1;
nt    = 216;
mut   = mu.^[1:nt]';
mutf  = mu.^[2:nt+1]';
for i=1:nsim;
  shocks   = randn(nt,7);
  [Ztf,Zt] = dlsim(G, H, G, H, shocks, [M;1]);

  Xt       = Zt(:,1:2);
  St       = Zt(:,3:9);
  Dt       = ([Xt,St]-ones(nt,1)*M')*[A,B]'+ones(nt,1)*log([hc,ic]);

  Xtf      = Ztf(:,1:2);
  Stf      = Ztf(:,3:9);
  Dtf      = ([Xtf,Stf]-ones(nt,1)*M')*[A,B]'+ones(nt,1)*log([hc,ic]);

  zt       = exp(St(:,1)).*mut;
  ztf      = exp(Stf(:,1)).*mutf;
  kct      = exp(Xt(:,1)).*mut;   
  kctf     = exp(Xtf(:,1)).*mutf;   
  knt      = exp(Xt(:,2)).*mut;   
  kntf     = exp(Xtf(:,2)).*mutf;   
  hct      = exp(Dt(:,1));       
  hctf     = exp(Dtf(:,1));       
  ict      = exp(Dt(:,2)).*mut;
  ictf     = exp(Dtf(:,2)).*mutf;
  tauct    = St(:,2);
  tauctf   = Stf(:,2);
  taunt    = St(:,3);
  tauntf   = Stf(:,3);
  tauht    = St(:,4);
  tauhtf   = Stf(:,4);
  gt       = exp(St(:,5)).*mut;
  gtf      = exp(Stf(:,5)).*mutf;
  thetat   = St(:,6);
  thetatf  = Stf(:,6);
  phit     = St(:,7);
  phitf    = Stf(:,7);

  hnt      = ( a2/a1* (1-phit)./(1-thetat).*((knt./zt).^(phit*rho)./ ...
               ((kct./zt).^(thetat*rho))).*hct.^(1-rho*(1-thetat)) ).^ ...
               (ones(nt,1)./(1-rho*(1-phit)));
  hntf     = ( a2/a1* (1-phitf)./(1-thetatf).*((kntf./ztf).^(phitf*rho)./ ...
               ((kctf./ztf).^(thetatf*rho))).*hctf.^(1-rho*(1-thetatf)) ).^ ...
               (ones(nt,1)./(1-rho*(1-phitf)));

  yct      = kct.^thetat.*(zt.*hct).^(1-thetat);
  yctf     = kctf.^thetatf.*(ztf.*hctf).^(1-thetatf);
  ynt      = knt.^phit  .*(zt.*hnt).^(1-phit);
  yntf     = kntf.^phitf.*(ztf.*hntf).^(1-phitf);
  yt       = (a1*yct.^rho+a2*ynt.^rho).^(1/rho);
  ytf      = (a1*yctf.^rho+a2*yntf.^rho).^(1/rho);
  pct      = a1 *(yt./yct).^(1-rho);
  pctf     = a1 *(ytf./yctf).^(1-rho);
  pnt      = a2 *(yt./ynt).^(1-rho);
  pntf     = a2 *(ytf./yntf).^(1-rho);
  rct      = pct.*thetat.*yct./kct;
  rctf     = pctf.*thetatf.*yctf./kctf;
  rnt      = pnt.*phit.*ynt./knt;
  rntf     = pntf.*phitf.*yntf./kntf;
  wt       = pct.*(1-thetat).*yct./hct;
  wtf      = pctf.*(1-thetatf).*yctf./hctf;
  ct       = (1-hct-hnt).*(1-tauht).*pct.*(1-thetat).*(kct./hct).^ ...
             thetat/psi .*zt.^(1-thetat);
  ctf      = (1-hctf-hntf).*(1-tauhtf).*pctf.*(1-thetatf).*(kctf./hctf).^ ...
             thetatf/psi .*ztf.^(1-thetatf);
  it       = yt-ct-gt;
  itf      = ytf-ctf-gtf;
  int      = it-ict;
  intf     = itf-ictf;
  dt       = pct.*yct-wt.*hct-ict;
  dtf      = pctf.*yctf-wtf.*hctf-ictf;
  Vct      = (rct+1-delc).*kct;
  Vctf     = (rctf+1-delc).*kctf;
  Vctx     = kctf;
  Vctxf    = [kctf(2:nt);kctf(nt)*mu];
  Ret      = (dt+Vctf)./Vct-1;
  Retx     = (dtf+Vctxf)./Vctx-1;
  Rete     = rctf-delc;

  Re(i)    = mean(Rete);

end;

[Erf,mean(Ret)-Erf,mean(Retx)-Erf,mean(Rete)-Erf]*400

mean((Re-Erf)*400)
