%GMOD_UC    computes equilibria for a growth model with two sectors --
%           a corporate sector and a noncorporate sector.  This code
%           allows for:
%              (1) adjustment costs on capital 
%              (2) unmeasured capital
%              (3) population growth.
%
%           Households solve:
%
%               max  E [sum  beta^t U(c[t],l[t]) N[t]
%              c,l,x
%
%             s.t. (1+tau_c[t]) c[t]+x_m1[t]+x_u1[t]+x_2[t] 
%                     = r_m1[t]k_m1[t]+r_u1[t]k_u1[t]+r_2[t]k_2[t]+w[t]n[t] 
%                           -tau_1[t]{(r_m1[t]-delm1-tau_k1[t])k_m1[t]
%                                     +r_u1[t]k_u1[t]-x_u1[t]} 
%                           -tau_2[t]{r_2[t]-del2-tau_k2[t]} k_2[t] 
%                           -tau_k1[t]k_m1[t]-tau_k2[t]k_2[t]
%                           -tau_n[t] w[t]n[t] + T[t]
%
%                        k_m1[t+1] = [(1-delm1) k_m1[t] + x_m1[t]]/(1+eta)
%                        k_u1[t+1] = [(1-delu1) k_u1[t] + x_u1[t]]/(1+eta)
%                        k_2[t+1]  = [(1-del2) k_2[t] + x_2[t]]/(1+eta)
%                        n[t]+l[t] = 1 
%                        N[t+1]    = (1+eta) N[t]
%
%           Corporate firms solve:
%
%               max   p_1 y_1 -r_m1 k_m1 -r_u1 k_u1 - w n_1    
%             k_1,n_1     s.t. y=k_m1^thetm k_u1^thetu (z n_1)^(1-thetm-thetu)
%
%
%           Noncorporate firms solve:
%
%               max   p_2 y_2 -r_2 k_2 - w n_2
%             k_2,n_2     s.t. y=k_2^phi (z n_2)^(1-phi)
%
%           Additionally,
%
%               y[t] = c[t] + x_m1[t] + x_u1[t] + x_2[t] + g[t]
%               n[t] = n_1[t] + n_2[t]
%               g[t] + T[t] = tau_1[t] (r_m1[t]-delm1) k_m1[t]
%                            +tau_2[t] (r_2[t]-del2) k_2[t] + tau_n[t]w[t]n[t]
%
%           Stochastic variables:  
%               ln z[t+1]  = ln z[t] + ln (1+gam) + omega[t+1]
%               omega[t+1] = rho_o omega[t] + eps_o[t+1]
%               tau_c[t+1] = (1-rho_c) mtauc  + rho_c tau_c[t] + eps_c[t+1]
%               tau_1[t+1] = (1-rho_1) mtau1  + rho_1 tau_1[t] + eps_1[t+1]
%               tau_2[t+1] = (1-rho_2) mtau2  + rho_2 tau_2[t] + eps_2[t+1]
%               tau_n[t+1] = (1-rho_n) mtaun  + rho_n tau_n[t] + eps_n[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]
%               tau_k1[t+1]= (1-rho_x) mtauk1 + rho_x tau_k1[t]+ eps_x[t+1]
%               tau_k2[t+1]= (1-rho_y) mtauk2 + rho_y tau_k2[t]+ eps_y[t+1]
%
%           Computation: find x_m1(k_m1,k_u1,k_2,s), x_u1(k_m1,k_u1,k_2,s),
%                        and n_1(k_m1,k_u1,k_2,s) where
%                        s = [omega,tau_c,tau_1,tau_2,tau_n,ln g, ...
%                             thetm,thetu,phi,tau_k1,tau_k2]'
%                        NOTE: all nonstationary variables must be divided
%                              by z[t] before computing. 
%             
%

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

%______________________________________________________________________________
%
%PARAMETERS:
%
%  Preferences:
%
%    U(c,l)   = (c l^psi)^(1-sigma)/(1-sigma)
%    beta     = the discount factor 
%   
%  Technology:
%
%    y_1      = k_m1^thetam *k_u1^thetau * (z n_1)^(1-thetam-thetau)  
%                                                   = corporate production
%    y_2      = k_2^phi   * (z n_2)^(1-phi)         = noncorporate production
%    y        = A*(a1 y_1^rho + a2 y_2^rho)^(1/rho) = GNP
%
%  Capital accumulation:
%
%    k_m1[t+1] = [(1-delm1) k_m1[t] + x_m1[t]]/(1+eta)
%    k_u1[t+1] = [(1-delu1) k_u1[t] + x_u1[t]]/(1+eta)
%    k_2[t+1]  = [(1-del2)  k_2[t]  + x_2[t]]/(1+eta)
%
%  Other:
%
%    gamma     = growth rate of technology
%    eta       = growth rate of population
%
%______________________________________________________________________________
%

nonlin   = input('Calculate hours nonlinearly? (1=yes,0=no) ');

% Free parameters:

b        = .12;
bm1      = b;
bu1      = b;
b2       = b;
rho      = -2;
sigma    = 1.5;
mtau2    = .0;
mtaun    = .25;
delu1    = .0;
Mxu1     = .019;           % unmeasured investment

rhoo     = .0;
rhoc     = .0;
rho1     = .0;
rho2     = .0;
rhon     = .0;
rhog     = .0;
rhom     = .0;
rhou     = .0;
rhop     = .0;
rhox     = .0;
rhoy     = .0;

sigo     = 0.0133;
sigc     = 0.0;
sig1     = 0.0;
sig2     = 0.0;
sign     = 0.0;
corno    = 0.0;
sigg     = 0.0;
sigm     = 0.0;
sigu     = 0.0;
sigp     = 0.0;
sigx     = 0.0;
sigy     = 0.0;

P        = [rhoo,0,0,0,0,0,0,0,0,0,0;
            0,rhoc,0,0,0,0,0,0,0,0,0;
            0,0,rho1,0,0,0,0,0,0,0,0;
            0,0,0,rho2,0,0,0,0,0,0,0;
            0,0,0,0,rhon,0,0,0,0,0,0;
            0,0,0,0,0,rhog,0,0,0,0,0;
            0,0,0,0,0,0,rhom,0,0,0,0;
            0,0,0,0,0,0,0,rhou,0,0,0;
            0,0,0,0,0,0,0,0,rhop,0,0;
            0,0,0,0,0,0,0,0,0,rhox,0;
            0,0,0,0,0,0,0,0,0,0,rhoy];
Sige     = [sigo,0,0,0,corno*sigo/sqrt(1-corno*corno),0,0,0,0,0,0;
            0,sigc,0,0,0,0,0,0,0,0,0;
            0,0,sig1,0,0,0,0,0,0,0,0;
            0,0,0,sig2,0,0,0,0,0,0,0;
            0,0,0,0,sign/sqrt(1-corno*corno),0,0,0,0,0,0;
            0,0,0,0,0,sigg,0,0,0,0,0;
            0,0,0,0,0,0,sigm,0,0,0,0;
            0,0,0,0,0,0,0,sigu,0,0,0;
            0,0,0,0,0,0,0,0,sigp,0,0;
            0,0,0,0,0,0,0,0,0,sigx,0;
            0,0,0,0,0,0,0,0,0,0,sigy];

% From data 

Dn       = .25;                % assume working 25% of time
Dgamma   = .00495;             % qtrly growth in technology
Deta     = .00245;             % qtrly growth in population
Dc       = .544;               % private consumption
Dxm1     = .100;               % corporate tangible investment
Dx2      = .207;               % noncorporate tangible investment
Dg       = .180;               % government consumption 

Dw1      = .378;               % corporate employee compensation
Dw2      = .246;               % noncorporate employee compensation
Dr1x     = .162;               % VA corp-W corp-xu1
Dr2      = .245;               % VA noncorp-W noncorp

Dkm1     = 1.042*4;            % corporate tangible assets
Dk2      = 2.447*4;            % noncorporate tangible assets

mtauc    = .047/Dc;            % consumption tax rate
mtau1    = .026/.073;          % corporate profit tax rate
mtauk1   = .02/Dkm1;           % corporate property tax rate 
mtauk2   = .012/Dk2;           % noncorporate property tax rate 



% Model parameters:

gamma    = Dgamma;
eta      = Deta;
Dku1     = Mxu1/(delu1+gamma+eta+gamma*eta);
delm1    = Dxm1/Dkm1 - (gamma+eta+gamma*eta);
del2     = Dx2/Dk2 - (gamma+eta+gamma*eta);
mphi     = Dr2/(Dw2+Dr2);
beta     = (1+gamma)^sigma/((1-mtau2)*(Dr2/Dk2-del2-mtauk2)+1);
betg     = beta*(1+gamma)^(-sigma);
Mru1     = (1-mtau2)*(Dr2/Dk2-del2-mtauk2)+delu1;
Mrm1     = (Mru1-delu1)/(1-mtau1)+delm1+mtauk1;
Dy       = Dc+Dxm1+Mxu1+Dx2+Dg;
Dva1     = Dr1x+Mxu1+Dw1;
Dva2     = Dr2+Dw2;
mthetm   = Mrm1*Dkm1/Dva1;
mthetu   = Mru1*Dku1/Dva1;
Mw       = ((1-mphi)*Dva2+(1-mthetm-mthetu)*Dva1)/Dn;
Mn1      = Dw1/Mw;
Mn2      = Dw2/Mw;
psi      = (1-mtaun)*Mw*(1-Mn1-Mn2)/(1+mtauc)/Dc;
My1      = Dkm1^mthetm*Dku1^mthetu*Mn1^(1-mthetm-mthetu);
My2      = Dk2^mphi*Mn2^(1-mphi);
Mp1      = Dva1/My1;
Mp2      = Dva2/My2;
Aa1      = Mp1*(Dy/My1)^(rho-1);
Aa2      = Mp2*(Dy/My2)^(rho-1);
scale    = (Aa1+Aa2)^(1/rho);
a1       = Aa1/scale^rho;
a2       = Aa2/scale^rho;
grate    = gamma+eta+eta*gamma;
mg       = Dg;

param    = [a1;a2;bm1;bu1;b2;beta;delm1;delu1;del2;eta;gamma;mphi;
            mthetm;mthetu;psi;rho;scale;sigma];
mlg      = log(mg);

%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

y1y    = My1/Dy;  % Guess of corporate GNP to total GNP
n1     = Mn1;     % Guess of hours of work in the corporate sector

z0     = [y1y;n1];
del    = diag(z0) * .00001;

for i=1:20;

   y1y = z0(1); 
   n1  = z0(2); 
   y2y = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
   p1  = scale^rho*a1*y1y^(rho-1);
   p2  = scale^rho*a2*y2y^(rho-1);
   r2  = (1/betg-1)/(1-mtau2)+del2+mtauk2;
   ru1 = 1/betg-1+delu1;
   rm1 = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
   nkm1= (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
         (1/(1-mthetm-mthetu)); 
   nk2 = (r2/mphi/p2)^(1/(1-mphi));
   km1 = n1/nkm1;
   ku1 = mthetu*rm1*km1/mthetm/ru1;
   nku1= n1/ku1;
   y1  = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
   y   = y1/y1y;
   y2  = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2; 
   xm1 = (grate+delm1)*km1;
   xu1 = (grate+delu1)*ku1;
   x2  = (grate+del2)*k2;
   c   = y-xm1-xu1-x2-mg;

   f   = [p1*(1-mthetm-mthetu)*nkm1^(-mthetm)*nku1^(-mthetu)- ...
             p2*(1-mphi)*nk2^(-mphi);
          psi*c-(1-mtaun)*p2*(1-mphi)*nk2^(-mphi)*(1-n1-n2)/(1+mtauc)];

   y1y = z0(1) + del(1,1);
   y2y = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
   p1  = scale^rho*a1*y1y^(rho-1);
   p2  = scale^rho*a2*y2y^(rho-1);
   r2  = (1/betg-1)/(1-mtau2)+del2+mtauk2;
   ru1 = 1/betg-1+delu1;
   rm1 = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
   nkm1= (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
         (1/(1-mthetm-mthetu)); 
   nk2 = (r2/mphi/p2)^(1/(1-mphi));
   km1 = n1/nkm1;
   ku1 = mthetu*rm1*km1/mthetm/ru1;
   nku1= n1/ku1;
   y1  = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
   y   = y1/y1y;
   y2  = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2;
   xm1 = (grate+delm1)*km1;
   xu1 = (grate+delu1)*ku1;
   x2  = (grate+del2)*k2;
   c   = y-xm1-xu1-x2-mg;

   df(:,1) ...
       = ([p1*(1-mthetm-mthetu)*nkm1^(-mthetm)*nku1^(-mthetu)- ...
             p2*(1-mphi)*nk2^(-mphi);
           psi*c-(1-mtaun)*p2*(1-mphi)*nk2^(-mphi)*(1-n1-n2)/(1+mtauc)]- ...
           f)/del(1,1);

   y1y = z0(1);
   n1  = z0(2) + del(2,2);
   y2y = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
   p1  = scale^rho*a1*y1y^(rho-1);
   p2  = scale^rho*a2*y2y^(rho-1);
   r2  = (1/betg-1)/(1-mtau2)+del2+mtauk2;
   ru1 = 1/betg-1+delu1;
   rm1 = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
   nkm1= (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
         (1/(1-mthetm-mthetu)); 
   nk2 = (r2/mphi/p2)^(1/(1-mphi));
   km1 = n1/nkm1;
   ku1 = mthetu*rm1*km1/mthetm/ru1;
   nku1= n1/ku1;
   y1  = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
   y   = y1/y1y;
   y2  = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2;
   xm1 = (grate+delm1)*km1;
   xu1 = (grate+delu1)*ku1;
   x2  = (grate+del2)*k2;
   c   = y-xm1-xu1-x2-mg;

   df(:,2) ...
       = ([p1*(1-mthetm-mthetu)*nkm1^(-mthetm)*nku1^(-mthetu)- ...
             p2*(1-mphi)*nk2^(-mphi);
           psi*c-(1-mtaun)*p2*(1-mphi)*nk2^(-mphi)*(1-n1-n2)/(1+mtauc)]- ...
           f)/del(2,2);

   z0  = z0 - df\f;

end;
y1y  = z0(1);
n1   = z0(2);
y2y  = ((scale^(-rho)-a1*y1y^rho)/a2)^(1/rho);
p1   = scale^rho*a1*y1y^(rho-1);
p2   = scale^rho*a2*y2y^(rho-1);
r2   = (1/betg-1)/(1-mtau2)+del2+mtauk2;
ru1  = 1/betg-1+delu1;
rm1  = (1/betg-1)/(1-mtau1)+delm1+mtauk1;
nkm1 = (rm1/mthetm/p1 * (mthetm*ru1/mthetu/rm1)^mthetu )^ ...
       (1/(1-mthetm-mthetu)); 
nk2  = (r2/mphi/p2)^(1/(1-mphi));

km1  = n1/nkm1;
ku1  = mthetu*rm1*km1/mthetm/ru1;
y1   = km1^mthetm*ku1^mthetu*n1^(1-mthetm-mthetu);
y    = y1/y1y;
y2   = (y^rho/(scale^rho*a2)-a1/a2*y1^rho)^(1/rho);
k2   = y2*nk2^(mphi-1);
km1f = (1+gamma)*km1;
ku1f = (1+gamma)*ku1;
k2f  = (1+gamma)*k2;
n2   = nk2*k2;
xm1  = (grate+delm1)*km1;
xu1  = (grate+delu1)*ku1;
x2   = (grate+del2)*k2;
c    = y-xm1-xu1-x2-mg;

rf   = (1+gamma)^sigma/beta-1;
rm1  = p1*mthetm*y1/km1;
ru1  = p1*mthetu*y1/ku1;
r2   = p2*mphi*y2/k2;
w    = p2*(1-mphi)*y2/n2;  

va1  = p1*y1;
va2  = p2*y2;
e2   = (r2-del2-mtauk2)*k2;
ptax = mtauk1*km1;
ctax = mtau1*((rm1-delm1-mtauk1)*km1+ru1*ku1-xu1);
re   = xm1+xu1-delm1*km1;
div  = va1-w*n1-ptax-ctax-xu1-xm1;
e1   = va1-w*n1-delm1*km1-ptax;
ate1 = va1-w*n1-delm1*km1-ptax-ctax;
atax = mtauc*c+mtaun*(n1+n2)+ctax+ptax+mtau2*e2+mtauk2*k2;
y    = c+xm1+xu1+x2+mg;
ir   = rf*4;
dv   = km1f+(1-mtau1)*ku1f;
tv   = dv + .0156/rf;
dev  = (ate1-xu1)/dv;
tev  = (ate1-xu1+.0156)/tv;
gnp  = y+mtauc*c-xu1+.037-.024-.012-.063-.016;
ibt  = mtauc*c+mtauk1*km1+mtauk2*k2;


%______________________________________________________________________________



%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

% The ordering of variables to the residuals is
%
%     Z  =  [ k_{m1,t+1},   k_{u1,t+1},   k_{2,t+1}, 
%             x_{m1,t+1},   x_{u1,t+1},   n_{1,t+1},
%             k_{m1,t},     k_{u1,t},     k_{2,t}, 
%             x_{m1,t},     x_{u1,t},     n_{1,t},
%             omega_{t+1},  tau_{c,t+1},  tau_{1,t+1}, 
%             tau_{2,t+1},  tau_{h,t+1},  g_{t+1}, 
%             thetm_{t+1},  thetu_{t+1},  phi_{t+1}, 
%             tauk_{1,t+1}, tauk_{2,t+1}, 
%             omega_t,      tau_{c,t},    tau_{1,t}, 
%             tau_{2,t},    tau_{h,t},    g_t, 
%             thetm_t,      thetu_t,      phi_t
%             tauk_{1,t},   tauk_{2,t}]' 
%

Z        = [log([km1,ku1,k2,xm1,xu1,n1,km1,ku1,k2,xm1,xu1,n1]), ...
            0,mtauc,mtau1,mtau2,mtaun,mlg,mthetm,mthetu,mphi,mtauk1,mtauk2, ...
            0,mtauc,mtau1,mtau2,mtaun,mlg,mthetm,mthetu,mphi,mtauk1,mtauk2]';
nZ       = 34;
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_gmod_uc(Zp,param);
  [Rm,R2m] = res_gmod_uc(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_{m1,t} ]   D_t = [ x_{m1,t} ]   S_t = [ omega_t   ]
%              [ k_{u1,t} ]         [ x_{u1,t} ]         [ tau_{c,t} ]
%              [ k_{2,t}  ]         [ n_{1,t}  ]         [ tau_{1,t} ]
%                                                        [ tau_{2,t} ]
%                                                        [ tau_{n,t} ]
%                                                        [ g_t       ]
%                                                        [ thetm_t   ]
%                                                        [ thetu_t   ]
%                                                        [ phi_t     ]
%                                                        [ tauk_{1,t}]
%                                                        [ tauk_{2,t}]
%
%
% The solution to the deterministic part of the problem is given by
%
%     D_t = A X_t
%
%

Psi1  = dR(1:3,1:3);
Psi2  = dR(1:3,4:6);
Psi3  = dR(1:3,7:9);
Psi4  = dR(1:3,10:12);
Psi5  = dR(1:3,13:23);
Psi6  = dR(1:3,24:34);

Phi1  = dR(4:6,1:3);
Phi2  = dR(4:6,4:6);
Phi3  = dR(4:6,7:9);
Phi4  = dR(4:6,10:12);
Phi5  = dR(4:6,13:23);
Phi6  = dR(4:6,24:34);

alp1  = dR2(7:9);
alp2  = dR2(10:12);
alp3  = dR2(24:34);

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


%______________________________________________________________________________


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

vB    = (kron(P',Psi2)+kron(eye(11),-(Psi1+Psi2*A)*(Phi1\Phi4))+ ...
          kron(eye(11),Psi4))\vec((Psi1+Psi2*A)*(Phi1\Phi5)*P-Psi5*P+ ...
                                 (Psi1+Psi2*A)*(Phi1\Phi6)-Psi6);

B     = zeros(3,11);
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] ]
%______________________________________________________________________________


%______________________________________________________________________________
% 
% SIMULATIONS FOR THE 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('Simulating time series ...')

G     = [-Phi1\(Phi3+Phi4*A),-Phi1\(Phi5*P+Phi6+Phi4*B);
          zeros(11,3), P];
H     = [ -(Phi1\Phi5)* Sige;
          Sige
          zeros(1,11)];
M     = [log(km1);log(ku1);log(k2);0;mtauc;mtau1;mtau2;
         mtaun;mlg;mthetm;mthetu;mphi;mtauk1;mtauk2];
G     = [G, (eye(14)-G)*M;
         zeros(1,14),1];


alp   = [alp1+alp2*A,alp2*B+alp3]';
e     = zeros(14,1);
e(4)  = 1;
gam   = (sigma*rhoo*e'+alp'*(eye(14)-G(1:14,1:14)))';
kap   = (-sigma*e'+alp')';
Vx    = dlyap(G(1:14,1:14),H(1:14,:)*H(1:14,:)');
Erf   = exp((gam'*Vx*gam-kap'*H(1:14,:)*H(1:14,:)'*kap)/2)* ...
        (1+gamma)^sigma/beta-1;
Erf   = ((1+Erf)^4-1)*100;

randn('seed',10462)
nsim  = 100;
nt    = 212;   % make this a multiple of 4
nyr   = nt/4; 
S     = zeros(nsim,10);
C     = zeros(nsim,7);
AC    = zeros(nsim,5);
mRetx = zeros(nsim,1);
mVal  = zeros(nsim,1);

for i=1:nsim;
  shocks   = randn(nt+2,11);
  [Ztf,Zt] = dlsim(G, H, G, H, shocks, [M;1]);

  Xt       = Zt(:,1:3);
  St       = Zt(:,4:14);
  Dt       = ([Xt,St]-ones(nt+2,1)*M')*[A,B]'+ones(nt+2,1)*log([xm1,xu1,n1]);
  omegat   = St(:,1);
  zt       = dlsim([1,log(1+gamma);0,1],[1;0],[1,log(1+gamma)],1,omegat,[0;1]);
  zt       = exp(zt)/(1+gamma);

  km1t     = exp(Xt(:,1)).*zt;
  km1tf    = km1t(2:nt+1);
  km1tf2   = km1t(3:nt+2);
  km1t     = km1t(1:nt);
  ku1t     = exp(Xt(:,2)).*zt;
  ku1tf    = ku1t(2:nt+1);
  ku1tf2   = ku1t(3:nt+2);
  ku1t     = ku1t(1:nt);
  k2t      = exp(Xt(:,3)).*zt;
  k2tf     = k2t(2:nt+1);
  k2tf2    = k2t(3:nt+2);
  k2t      = k2t(1:nt);
  ztf      = zt(2:nt+1);
  zt       = zt(1:nt);
  dhat     = delm1+grate;
  xm1t     = (1+bm1*dhat-sqrt(1+2*bm1*(1+grate)-2*bm1*(1+eta)*km1tf./ ...
                km1t))/bm1 .* km1t;
  xm1tf    = (1+bm1*dhat-sqrt(1+2*bm1*(1+grate)-2*bm1*(1+eta)*km1tf2./ ...
                km1tf))/bm1 .* km1tf;
  dhat     = delu1+grate;
  xu1t     = (1+bu1*dhat-sqrt(1+2*bu1*(1+grate)-2*bu1*(1+eta)*ku1tf./ ...
                ku1t))/bu1 .* ku1t;
  xu1tf    = (1+bu1*dhat-sqrt(1+2*bu1*(1+grate)-2*bu1*(1+eta)*ku1tf2./ ...
                ku1tf))/bu1 .* ku1tf;
  dhat     = del2+grate;
  x2t      = (1+b2*dhat-sqrt(1+2*b2*(1+grate)-2*b2*(1+eta)*k2tf./k2t))/ ...
                b2 .* k2t;
  x2tf     = (1+b2*dhat-sqrt(1+2*b2*(1+grate)-2*b2*(1+eta)*k2tf2./k2tf))/ ...
                b2 .* k2tf;
  xt       = xm1t+xu1t+x2t;
  xtf      = xm1tf+xu1tf+x2tf;
  tauct    = St(1:nt,2);
  tauctf   = St(2:nt+1,2);
  tau1t    = St(1:nt,3);
  tau1tf   = St(2:nt+1,3);
  tau2t    = St(1:nt,4);
  tau2tf   = St(2:nt+1,4);
  taunt    = St(1:nt,5);
  tauntf   = St(2:nt+1,5);
  gt       = exp(St(1:nt,6)).*zt;
  gtf      = exp(St(2:nt+1,6)).*ztf;
  thetmt   = St(1:nt,7);
  thetmtf  = St(2:nt+1,7);
  thetut   = St(1:nt,8);
  thetutf  = St(2:nt+1,8);
  phit     = St(1:nt,9);
  phitf    = St(2:nt+1,9);
  tauk1t   = St(1:nt,10);
  tauk1tf  = St(2:nt+1,10);
  tauk2t   = St(1:nt,11);
  tauk2tf  = St(2:nt+1,11);

  n1t      = exp(Dt(1:nt,3));
  n1tf     = exp(Dt(2:nt+1,3));
  if nonlin;
    n1t      = secant('foc_uc',n1t,[km1t,ku1t,k2t,zt,tauct,taunt,gt, ...
                       thetmt,thetut,phit,xt, ...
                       ones(nt,1)*[a1,a2,rho,psi,scale]]);
    tem      = secant('foc_uc',exp(Dt(nt+1,3)),[km1tf(nt),ku1tf(nt), ...
                       k2tf(nt),ztf(nt),tauctf(nt),tauntf(nt), ...
                       gtf(nt),thetmtf(nt),thetutf(nt),phitf(nt), ...
                       xtf(nt),a1,a2,rho,psi,scale]);
    n1tf     = [n1t(2:nt);tem];
  end;

  n2t      = ( a2/a1* (1-phit)./(1-thetmt-thetut).*((k2t./zt).^(phit*rho)./ ...
               ( (km1t./zt).^(thetmt*rho).*(ku1t./zt).^(thetut*rho))).* ...
                  n1t.^(1-rho*(1-thetmt-thetut)) ).^(ones(nt,1)./ ...
                 (1-rho*(1-phit)));
  n2tf     = ( a2/a1* (1-phitf)./(1-thetmtf-thetutf).*((k2tf./ztf).^ ...
                 (phitf*rho)./( (km1tf./ztf).^(thetmtf*rho).* ...
                 (ku1tf./ztf).^(thetutf*rho))).*n1tf.^ ...
                 (1-rho*(1-thetmtf-thetutf)) ).^(ones(nt,1)./(1-rho*(1-phitf)));
  y1t      = km1t.^thetmt.*ku1t.^thetut.*(zt.*n1t).^(1-thetmt-thetut);
  y1tf     = km1tf.^thetmtf.*ku1tf.^thetutf.*(ztf.*n1tf).^(1-thetmtf-thetutf);
  y2t      = k2t.^phit  .*(zt.*n2t).^(1-phit);
  y2tf     = k2tf.^phitf.*(ztf.*n2tf).^(1-phitf);
  yt       = scale*(a1*y1t.^rho+a2*y2t.^rho).^(1/rho);
  ytf      = scale*(a1*y1tf.^rho+a2*y2tf.^rho).^(1/rho);
  ct       = yt-xt-gt;
  ctf      = ytf-xtf-gtf;
  p1t      = scale^rho*a1 *(yt./y1t).^(1-rho);
  p1tf     = scale^rho*a1 *(ytf./y1tf).^(1-rho);
  p2t      = scale^rho*a2 *(yt./y2t).^(1-rho);
  p2tf     = scale^rho*a2 *(ytf./y2tf).^(1-rho);
  rm1t     = p1t.*thetmt.*y1t./km1t;
  rm1tf    = p1tf.*thetmtf.*y1tf./km1tf;
  ru1t     = p1t.*thetut.*y1t./ku1t;
  ru1tf    = p1tf.*thetutf.*y1tf./ku1tf;
  r2t      = p2t.*phit.*y2t./k2t;
  r2tf     = p2tf.*phitf.*y2tf./k2tf;
  wt       = p2t.*(1-phit).*y2t./n2t; 
  wtf      = p2tf.*(1-phitf).*y2tf./n2tf; 
  e2t      = (r2t-del2-tauk2t).*k2t;
  e2tf     = (r2tf-del2-tauk2tf).*k2tf;
  ctaxt    = tau1t.*((rm1t-delm1-tauk1t).*km1t+ru1t.*ku1t-xu1t);
  ctaxtf   = tau1tf.*((rm1tf-delm1-tauk1tf).*km1tf+ru1tf.*ku1tf-xu1tf);
  divt     = p1t.*y1t-wt.*n1t-tauk1t.*km1t-ctaxt-xu1t-xm1t;
  divtf    = p1tf.*y1tf-wtf.*n1tf-tauk1tf.*km1tf-ctaxtf-xu1tf-xm1tf;
  ate1t    = p1t.*y1t-wt.*n1t-delm1*km1t-tauk1t.*km1t-ctaxt;
  ate1tf   = p1tf.*y1tf-wtf.*n1tf-delm1*km1tf-tauk1tf.*km1tf-ctaxtf;

  it       = ones(nt,1)*Erf/100;
  itf      = ones(nt,1)*Erf/100;
  gnpt     = yt+tauct.*ct-xu1t+.037*yt/y-.036*it.*k2t/(Erf*k2/100)-.079*k2t/k2;
  vt       = km1tf./(1-bm1*(xm1t./km1t-grate-delm1))+ ....
             (1-tau1t).*ku1tf./(1-bu1*(xu1t./ku1t-grate-delu1));
  vtf      = km1tf2./(1-bm1*(xm1tf./km1tf-grate-delm1))+ ....
             (1-tau1tf).*ku1tf2./(1-bu1*(xu1tf./ku1tf-grate-delu1));
  tvt      = vt+4*.0156*(ate1t-xu1t+ctaxt)./it/(ate1-xu1+ctax);
  Retx     = (divtf+(1+eta)*vtf)./vt-1;
  tem      = cumprod(1+Retx);
  mRetx(i) = (tem(length(tem))^(1/nyr)-1)*100;
  mVal(i)  = mean(tvt./gnpt/4);

  dyt      = log((1+tauct).*ct+xm1t+x2t+gt)- ...
             hptrend(log((1+tauct).*ct+xm1t+x2t+gt),1600);
  dct      = log((1+tauct).*ct)-hptrend(log((1+tauct).*ct),1600);
  dxt      = log(xm1t+x2t)-hptrend(log(xm1t+x2t),1600);
  dgt      = log(gt)-hptrend(log(gt),1600);
  dnt      = log(n1t+n2t)-hptrend(log(n1t+n2t),1600);
  dwtax    = log(taunt.*wt.*(n1t+n2t))-hptrend(log(taunt.*wt.*(n1t+n2t)),1600);
  dktax    = log(ctaxt)-hptrend(log(ctaxt),1600);
  dctax    = log(tauct.*ct)-hptrend(log(tauct.*ct),1600);
  date1t   = log(ate1t-xu1t)-hptrend(log(ate1t-xu1t),1600);
  dxm1t    = zeros(4,nyr); dx2t=dxm1t;
  dxm1t(:) = xm1t;  dxm1t=sum(dxm1t)';
  dx2t(:)  = x2t;   dx2t =sum(dx2t)';
  dxm1t    = log(dxm1t)-hptrend(log(dxm1t),100);
  dx2t     = log(dx2t)-hptrend(log(dx2t),100);
  S(i,:)   = [std(dyt),std(dct),std(dxt),std(dgt),std(dnt), ...
              std(dwtax),std(dktax),std(dctax),std(date1t),std(thetmt)]*100;
  tem1     = corrcoef([dyt(1:nt-1),dyt(2:nt)]);
  tem2     = corrcoef([dct(1:nt-1),dct(2:nt)]);
  tem3     = corrcoef([dxt(1:nt-1),dxt(2:nt)]);
  tem4     = corrcoef([dgt(1:nt-1),dgt(2:nt)]);
  tem5     = corrcoef([date1t(1:nt-1),date1t(2:nt)]);
  AC(i,:)  = [tem1(1,2),tem2(1,2),tem3(1,2),tem4(1,2),tem5(1,2)];
  tem1     = corrcoef([thetmt(1:nt-1),thetmt(2:nt)]);
  tem2     = corrcoef([thetmt,dyt]);
  tem3     = corrcoef([date1t,dct]);
  tem4     = corrcoef([dxm1t,dx2t]);
  tem5     = corrcoef([dwtax,dyt]);
  tem6     = corrcoef([dktax,dyt]);
  tem7     = corrcoef([dctax,dyt]);
  C(i,:)   = [tem1(1,2),tem2(1,2),tem3(1,2),tem4(1,2), ...
              tem5(1,2),tem6(1,2),tem7(1,2)];

  disp(sprintf(['Simulation %g complete, Equity return = %5.2f,', ...
                ' RF rate = %5.2f, Val = %5.2f'],[i,mRetx(i),Erf,mVal(i)]))
end;

disp(' ')
disp(' ')
disp(' ')
disp('MEANS OF RATIOS (NOT DIVIDED BY GNP)')
tab1  = [gnp,y,c,xm1,x2,mg,xu1,mtauc*c,gnp-y-mtauc*c+xu1, ...
          gnp,delm1*km1+del2*k2,gnp-delm1*km1-del2*k2, ...
          ibt,gnp-delm1*km1-del2*k2-ibt, ...
          y,delm1*km1+del2*k2,mtauk1*km1+mtauk2*k2,w*(n1+n2),e2,e1, ...
          va1,w*n1,delm1*km1,ptax,e1,ctax,e1-ctax,div,re, ...
          va2,w*n2,del2*k2,mtauk2*k2,e2, ...
          (km1f+ku1f+k2f)/4,km1f/4,ku1f/4,k2f/4, ...
          xm1+xu1+x2,xm1,xu1,x2, ...
          delm1*km1+delu1*ku1+del2*k2,delm1*km1,delu1*ku1,del2*k2, ...
          atax,mtauc*c,mtaun*(n1+n2),ctax,ptax,mtau2*e2,mtauk2*k2, ...
          ([(1+gamma)^sigma/beta]^4-1)*100, ...
          dv/4,tv/4,dev*400,tev*400];

str1  = ['  Components of NIPA GNP, y+tau_c*c-x_u1-adjustments %6.3f\n';
         '    Model GNP:                                       %6.3f\n';
         '      Private Consumption less tax, c                %6.3f\n';
         '      Corporate measured investment, x_m1            %6.3f\n';
         '      Noncorporate investment, x_2                   %6.3f\n';
         '      Government Purchases, g                        %6.3f\n';
         '    PLUS: Consumption tax, tau_c*c                   %6.3f\n';
         '    LESS: Corporate unmeasured investment, xu1       %6.3f\n';
         '    LESS: Adjustments not made by NIPA               %6.3f\n';
         '                                                          \n';
         '  NIPA GNP --> Income                                %6.3f\n';
         '    - NIPA Depreciation, del_m1 *k_m1 + del_2 *k_2   %6.3f\n';
         '    = Net Product                                    %6.3f\n';
         '    - IBT, tau_c*c+tau_k1 k_m1+tau_k2 k_m2           %6.3f\n';
         '    = Income                                         %6.3f\n';
         '                                                          \n';
         '  Components of Model Income, p_1*y_1+p_2*y_2        %6.3f\n';
         '    Depreciation, delm1*km1+del2*k2                  %6.3f\n';
         '    Property tax, tau_k1*k_m1+tau_k2*k_2             %6.3f\n';
         '    Wages, w*(n_1+n_2)                               %6.3f\n';
         '    Noncorporate Earnings, (r_2-del_2-tau_k2)*k_2    %6.3f\n';
         '    Corporate Earnings (detail below)                %6.3f\n';
         '                                                          \n';
         '  Corporate Income, p_1*y_1                          %6.3f\n';
         '    Corporate wages, w*n_1                           %6.3f\n';
         '    Measured depreciation, del_m1*k_m1               %6.3f\n';
         '    Property tax, tau_k1*k_m1                        %6.3f\n';
         '    Corporate earnings                               %6.3f\n';
         '      Tax, tau_1*[(r_m1-del_m1-tau_k1)*k_m1               \n';
         '                  +r_u1*k_u1-x_u1]                   %6.3f\n';
         '      Earnings after Tax                             %6.3f\n';
         '        Dividends, p_1*y_1-wage-taxes-x_m1-x_u1      %6.3f\n';
         '        Retained earnings, x_m1+x_u1-del_m1*k_m1     %6.3f\n';
         '                                                          \n';
         '  Noncorporate Income, p_2*y_2                       %6.3f\n';
         '    Noncorporate wages, w*n_2                        %6.3f\n';
         '    Depreciation, del_2*k_2                          %6.3f\n';
         '    Property tax, tau_k2*k_2                         %6.3f\n';
         '    Noncorporate earnings, (r_2-del_2-tau_k2)*k_2    %6.3f\n';
         '                                                          \n';
         '  Capital, end of year                               %6.3f\n';
         '    Corporate - measured, k_m1''/4                    %6.3f\n';
         '    Corporate - unmeasured, k_u1''/4                  %6.3f\n';
         '    Noncorporate, k_2''/4                             %6.3f\n';
         '                                                          \n';
         '  Investment                                         %6.3f\n';
         '    Corporate - measured, x_m1                       %6.3f\n';
         '    Corporate - unmeasured, x_u1                     %6.3f\n';
         '    Noncorporate, x_2                                %6.3f\n';
         '                                                          \n';
         '  Depreciation                                       %6.3f\n';
         '    Corporate - measured, del_m1*k_m1                %6.3f\n';
         '    Corporate - unmeasured, del_u1*k_u1              %6.3f\n';
         '    Noncorporate, del_2*k_2                          %6.3f\n';
         '                                                          \n';
         '  Tax Revenues                                       %6.3f\n';
         '    Consumption taxes, tau_c*c                       %6.3f\n';
         '    Wage taxes, tau_n*w*(n_1+n_2)                    %6.3f\n';
         '    Corporate profit taxes, (see formula above)      %6.3f\n';
         '    Corporate property tax, tau_k1*k_m1              %6.3f\n';
         '    Noncorporate taxes, tau_2*e_2                    %6.3f\n';
         '    Noncorporate property tax, tau_k2*k_2            %6.3f\n';
         '                                                          \n';
         '  Interest rate, i=([(1+gamma)^sigma/beta]^4-1)*100  %6.3f\n';
         '                                                          \n';
         '  Value of Corporate Equities                             \n';
         '    Domestic, (k_m1''+(1-tau_1)k_u1'')/4               %6.3f\n';
         '    Total, (k_m1''+(1-tau_)k_u1'')/4+1.2/i             %6.3f\n';
         '                                                          \n';
         '  After-tax earnings/Corporate value                      \n';
         '    Domestic ATE/Domestic value                      %6.3f\n';
         '    Total ATE/Total value                            %6.3f\n'];

disp(sprintf(str1',tab1))
disp(' ')

disp('STANDARD DEVIATIONS OF LOGGED AND HP-FILTERED SERIES (%)')
str2   =['  Measured output, (1+tau_c)*c+x_m1+x_2+g            %6.3f\n';
         '  Consumption, (1+tau_c)*c                           %6.3f\n';
         '  Investment, x_m1+x_2                               %6.3f\n';
         '  Government Purchases, g                            %6.3f\n';
         '  Hours of work, n                                   %6.3f\n';
         '  Wage Taxes, tau_n*w*(n_1+n_2)                      %6.3f\n';
         '  Corporate taxes, tau_1*e_1                         %6.3f\n';
         '  Consumption taxes, tau_c*c                         %6.3f\n';
         '  After-tax earnings                                 %6.3f  '];
tab2   = [mean(S(:,1)),mean(S(:,2)),mean(S(:,3)),mean(S(:,4)), ...
          mean(S(:,5)),mean(S(:,6)),mean(S(:,7)),mean(S(:,8)),mean(S(:,9))];
disp(sprintf(str2',tab2))
disp(' ')

disp('AUTOCORRELATIONS OF LOGGED AND HP-FILTERED SERIES')
str2b  =['  Measured output, (1+tau_c)*c+x_m1+x_2+g            %6.3f\n';
         '  Consumption, (1+tau_c)*c                           %6.3f\n';
         '  Investment, x_m1+x_2                               %6.3f\n';
         '  Government Purchases, g                            %6.3f\n';
         '  After-tax earnings                                 %6.3f  '];
tab2b  = [mean(AC(:,1)),mean(AC(:,2)),mean(AC(:,3)),mean(AC(:,4)), ...
          mean(AC(:,5))];
disp(sprintf(str2b',tab2b))
disp(' ')

disp('CORRELATIONS FOR LOGGED AND HP-FILTERED VARIABLES')
str3   =['  Consumption & After-Tax earnings                   %6.3f\n';
         '  Corporate & Noncorporate Investment (annual)       %6.3f\n'; 
         '  GNP & Wage Taxes, tau_n*w*(n_1+n_2)                %6.3f\n';
         '  GNP & Corporate taxes, tau_1*e_1                   %6.3f\n';
         '  GNP & Consumption taxes, tau_c*c                   %6.3f\n'];
tab3   = [mean(C(:,3)),mean(C(:,4)),mean(C(:,5)),mean(C(:,6)),mean(C(:,7))];
disp(sprintf(str3',tab3))
disp(' ')

disp('STATISTICS FOR MEASURED CORPORATE CAPITAL SHARE')
str4  = ['  Mean                                               %6.3f\n';
         '  Standard Deviation (percent)                       %6.3f\n';
         '  Autocorrelation                                    %6.3f\n';
         '  Correlation with logged and filtered real GNP      %6.3f\n'];
tab4  = [mthetm,mean(S(:,10)),mean(C(:,1)),mean(C(:,2))];
disp(sprintf(str4',tab4))
disp(' ')

disp('MEAN RETURNS/VALUES')
str5  = ['  Risk-free rate                                     %6.3f\n';
         '  Corporate equity return                            %6.3f\n';
         '  Equity premium                                     %6.3f\n';
         '  Stock Value to GNP                                 %6.3f\n'];
tab5  = [Erf,mean(mRetx),mean(mRetx)-Erf,mean(mVal)];
disp(sprintf(str5',tab5))

