%GMOD_AC    computes equilibria for a growth model with two sectors --
%           a corporate sector and a noncorporate sector.  Adjustment
%           costs on capital are included.
%
%           Households solve:
%
%               max  E [sum  beta^t U(c[t],l[t])
%              c,l,x
%                   s.t. (1+tau_c[t]) c[t]+x_1[t]+x_2[t] 
%                            = (1-tau_1[t]) r_1[t]k_1[t]+tau_1[t] del1 k_1[t]
%                             +(1-tau_2[t]) r_2[t]k_2[t]+tau_2[t] del2 k_2[t]
%                             +(1-tau_n[t]) w[t] n[t] + T[t]
%                        k_1[t+1] = (1-del1) k_1[t] + x_1[t]
%                        k_2[t+1] = (1-del2) k_2[t] + x_2[t]
%                        n[t]+l[t]= 1 
%
%           Corporate firms solve:
%
%               max   p_1 y_1 -r_1 k_1 - w n_1    
%             k_1,n_1       s.t. y=k_1^theta (z n_1)^(1-theta)
%
%
%           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_1[t] + x_2[t] + g[t]
%               n[t] = n_1[t] + n_2[t]
%               g[t] + T[t] = tau_1[t] (r_1[t]-del1) k_1[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 mu + eta[t+1]
%               eta[t+1]   = rho_e eta[t] + eps_e[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]
%
%           Computation: find n_1(k_1,k_2,s), x_1(k_1,k_2,s)
%                        where s = [eta,tau_c,tau_1,tau_2,tau_n,ln g,theta,phi]'
%                        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_1^theta * (z n_1)^(1-theta)       = corporate production
%    y_2      = k_2^phi   * (z n_2)^(1-phi)         = noncorporate production
%    y        = (a1 y_1^rho + a2 y_2^rho)^(1/rho)   = GNP
%
%  Capital accumulation
%
%    k_1[t+1] = (1-del1) k_1[t] + x_1[t]
%    k_2[t+1] = (1-del2) k_2[t] + x_2[t]
%
%______________________________________________________________________________
%

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

psi      = .01;
sigma    = 5;
beta     = .99*(1.005)^(sigma-1); 
a1       = 1;
a2       = 1;
b1       = .00001;
b2       = .00001;
rho      = .999;
del1     = .025;
del2     = .025;
mu       = 1.005;

mtauc    = .0;
mtau1    = .0;
mtau2    = .0; 
mtaun    = .0;
mg       = .00001;
mthet    = .36;
mphi     = .36;

rhoe     = .0;
rhoc     = .0;
rho1     = .0;
rho2     = .0;
rhon     = .0;
rhog     = .0;
rhot     = .0;
rhop     = .0;

sige     = 0.018;
sigc     = 0.0;
sig1     = 0.0;
sig2     = 0.0;
sign     = 0.0;
sigg     = 0.0;
sigt     = 0.0;
sigp     = 0.0;
cortp    = 0.0;

param    = [a1;a2;b1;b2;beta;del1;del2;mphi;mthet;mu;psi;rho;sigma];
betm     = beta*mu^(-sigma);
mlg      = log(mg);

%______________________________________________________________________________



%______________________________________________________________________________
%
% STEADY STATE FOR THE ECONOMY
%

y1y    = .90;  % Guess of corporate GNP to total GNP
n1     = .22;  % 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 = ((1-a1*y1y^rho)/a2)^(1/rho);
   p1  = a1*y1y^(rho-1);
   p2  = a2*y2y^(rho-1);
   nk1 = ((1-betm+(1-mtau1)*betm*del1)/(betm*(1-mtau1)*mthet*p1))^(1/(1-mthet));
   nk2 = ((1-betm+(1-mtau2)*betm*del2)/(betm*(1-mtau2)*mphi*p2))^(1/(1-mphi));
   k1  = n1/nk1;
   y1  = nk1^(1-mthet)*k1;
   y   = y1/y1y;
   y2  = (y^rho/a2-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2; 
   x1  = (mu-1+del1)*k1;
   x2  = (mu-1+del2)*k2;
   c   = y-x1-x2-mg;
    
   f   = [p1*(1-mthet)*nk1^(-mthet)-p2*(1-mphi)*nk2^(-mphi);
          psi*c-(1-mtaun)*p1*(1-mthet)*nk1^(-mthet)*(1-n1-n2)/(1+mtauc)];

   y1y = z0(1) + del(1,1);
   y2y = ((1-a1*y1y^rho)/a2)^(1/rho);
   p1  = a1*y1y^(rho-1);
   p2  = a2*y2y^(rho-1);
   nk1 = ((1-betm+(1-mtau1)*betm*del1)/(betm*(1-mtau1)*mthet*p1))^(1/(1-mthet));
   nk2 = ((1-betm+(1-mtau2)*betm*del2)/(betm*(1-mtau2)*mphi*p2))^(1/(1-mphi));
   k1  = n1/nk1;
   y1  = nk1^(1-mthet)*k1;
   y   = y1/y1y;
   y2  = (y^rho/a2-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2; 
   x1  = (mu-1+del1)*k1;
   x2  = (mu-1+del2)*k2;
   c   = y-x1-x2-mg;
    
   df(:,1) ...
       = ([p1*(1-mthet)*nk1^(-mthet)-p2*(1-mphi)*nk2^(-mphi);
           psi*c-(1-mtaun)*p1*(1-mthet)*nk1^(-mthet)*(1-n1-n2)/(1+mtauc)]- ...
           f)/del(1,1);

   y1y = z0(1);
   n1  = z0(2) + del(2,2);
   y2y = ((1-a1*y1y^rho)/a2)^(1/rho);
   p1  = a1*y1y^(rho-1);
   p2  = a2*y2y^(rho-1);
   nk1 = ((1-betm+(1-mtau1)*betm*del1)/(betm*(1-mtau1)*mthet*p1))^(1/(1-mthet));
   nk2 = ((1-betm+(1-mtau2)*betm*del2)/(betm*(1-mtau2)*mphi*p2))^(1/(1-mphi));
   k1  = n1/nk1;
   y1  = nk1^(1-mthet)*k1;
   y   = y1/y1y;
   y2  = (y^rho/a2-a1/a2*y1^rho)^(1/rho);
   k2  = y2*nk2^(mphi-1);
   n2  = nk2*k2; 
   x1  = (mu-1+del1)*k1;
   x2  = (mu-1+del2)*k2;
   c   = y-x1-x2-mg;
    
   df(:,2) ...
       = ([p1*(1-mthet)*nk1^(-mthet)-p2*(1-mphi)*nk2^(-mphi);
           psi*c-(1-mtaun)*p1*(1-mthet)*nk1^(-mthet)*(1-n1-n2)/(1+mtauc)]- ...
           f)/del(2,2);

   z0  = z0 - df\f;

end;
y1y = z0(1);
n1  = z0(2);
y2y = ((1-a1*y1y^rho)/a2)^(1/rho);
p1  = a1*y1y^(rho-1);
p2  = a2*y2y^(rho-1);
nk1 = ((1-betm+(1-mtau1)*betm*del1)/(betm*(1-mtau1)*mthet*p1))^(1/(1-mthet));
nk2 = ((1-betm+(1-mtau2)*betm*del2)/(betm*(1-mtau2)*mphi*p2))^(1/(1-mphi));
k1  = n1/nk1;
y1  = nk1^(1-mthet)*k1;
y   = y1/y1y;
y2  = (y^rho/a2-a1/a2*y1^rho)^(1/rho);
k2  = y2*nk2^(mphi-1);
n2  = nk2*k2; 
x1  = (mu-1+del1)*k1;
x2  = (mu-1+del2)*k2;
x   = x1+x2;
k1f = mu*k1;
k2f = mu*k2;
c   = y-x1-x2-mg;
rf  = mu^sigma/beta-1;
r1  = p1*mthet*y1/k1;
r2  = p2*mphi*y2/k2;
w   = p1*(1-mthet)*y1/n1;  
e1  = (r1-del1)*k1;
e2  = (r2-del2)*k2;
rev = mtauc*c+mtau1*e1+mtau2*e2+mtaun*w*(n1+n2);
ate1= (1-mtau1)*e1;
d   = ate1-mu*k1+k1;
re  = mu*k1-k1;
v   = mu*k1;
dep1= del1*k1;
dep2= del2*k2;
dep = dep1+dep2;

%______________________________________________________________________________




%______________________________________________________________________________
%
% LINEARIZED FIRST ORDER CONDITIONS
%

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

Z        = [log([k1,k2,n1,x1,k1,k2,n1,x1]), ...
            0,mtauc,mtau1,mtau2,mtaun,mlg,mthet,mphi, ...
            0,mtauc,mtau1,mtau2,mtaun,mlg,mthet,mphi]';
nZ       = 24;
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_ac(Zp,param);
  [Rm,R2m] = res_gmod_ac(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_{1,t} ]   D_t = [ n_{1,t} ]   S_t = [ eta_t     ]
%              [ k_{2,t} ]         [ x_{1,t} ]         [ tau_{c,t} ]
%                                                      [ tau_{1,t} ]
%                                                      [ tau_{2,t} ]
%                                                      [ tau_{n,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:16);
Psi6  = dR(1:2,17:24);

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:16);
Phi6  = dR(3:4,17:24);

alp1  = dR2(5:6);
alp2  = dR2(7:8);
alp3  = dR2(17:24);

[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     = [rhoe,0,0,0,0,0,0,0;
         0,rhoc,0,0,0,0,0,0;
         0,0,rho1,0,0,0,0,0;
         0,0,0,rho2,0,0,0,0;
         0,0,0,0,rhon,0,0,0;
         0,0,0,0,0,rhog,0,0;
         0,0,0,0,0,0,rhot,0;
         0,0,0,0,0,0,0,rhop];
vB    = (kron(P',Psi2)+kron(eye(8),-(Psi1+Psi2*A)*(Phi1\Phi4))+ ...
          kron(eye(8),Psi4))\vec((Psi1+Psi2*A)*(Phi1\Phi5)*P-Psi5*P+ ...
                                 (Psi1+Psi2*A)*(Phi1\Phi6)-Psi6);

B     = zeros(2,8);
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('Simulating time series ...')

Sige  = [sige,0,0,0,0,0,0,0;
         0,sigc,0,0,0,0,0,0;
         0,0,sig1,0,0,0,0,0;
         0,0,0,sig2,0,0,0,0;
         0,0,0,0,sign,0,0,0;
         0,0,0,0,0,sigg,0,0;
         0,0,0,0,0,0,sigt,cortp*sigt/sqrt(1-cortp*cortp);
         0,0,0,0,0,0,0,sigp/sqrt(1-cortp*cortp)];
G     = [-Phi1\(Phi3+Phi4*A),-Phi1\(Phi5*P+Phi6+Phi4*B);
          zeros(8,2), P];
H     = [ -(Phi1\Phi5)* Sige;
          Sige
          zeros(1,8)];
M     = [log(k1);log(k2);0;mtauc;mtau1;mtau2;mtaun;mlg;mthet;mphi];
G     = [G, (eye(10)-G)*M;
         zeros(1,10),1];


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

randn('seed',10462)
nsim  = 1;
nt    = 212;   % make this a multiple of 4
nyr   = 53; 
S     = zeros(nsim,6);
C     = zeros(nsim,4);
mRetx = zeros(nsim,1);

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

  Xt       = Zt(:,1:2);
  St       = Zt(:,3:10);
  Dt       = ([Xt,St]-ones(nt+2,1)*M')*[A,B]'+ones(nt+2,1)*log([n1,x1]);
  etat     = St(:,1);
  zt       = dlsim([1,log(mu);0,1],[1;0],[1,log(mu)],1,etat,[0;1]);
  zt       = exp(zt);

  k1t      = exp(Xt(:,1)).*zt;
  k1tf     = k1t(2:nt+1);
  k1tf2    = k1t(3:nt+2);
  k1t      = k1t(1:nt);
  k2t      = exp(Xt(:,2)).*zt;
  k2tf     = k2t(2:nt+1);
  k2tf2    = k2t(3:nt+2);
  k2t      = k2t(1:nt);
  ztf      = zt(2:nt+1);
  zt       = zt(1:nt);
  x1t      = (1+b1*del1-sqrt(1+2*b1-2*b1*k1tf./k1t))/b1 .* k1t;
  x1tf     = (1+b1*del1-sqrt(1+2*b1-2*b1*k1tf2./k1tf))/b1 .* k1tf;
  x2t      = (1+b2*del2-sqrt(1+2*b2-2*b2*k2tf./k2t))/b2 .* k2t;
  x2tf     = (1+b2*del2-sqrt(1+2*b2-2*b2*k2tf2./k2tf))/b2 .* k2tf;
  xt       = x1t+x2t;
  xtf      = x1tf+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;
  thetat   = St(1:nt,7);
  thetatf  = St(2:nt+1,7);
  phit     = St(1:nt,8);
  phitf    = St(2:nt+1,8);

  n1t      = exp(Dt(1:nt,1));
  n1tf     = exp(Dt(2:nt+1,1));
  if nonlin;
    n1t      = secant('foc',n1t,[k1t,k2t,zt,tauct,tau1t,tau2t,taunt,gt, ...
                       thetat,phit,xt,ones(nt,1)*[a1,a2,rho,psi]]);
    tem      = secant('foc',exp(Dt(nt+1,1)),[k1tf(nt),k2tf(nt),ztf(nt), ...
                      tauctf(nt),tau1tf(nt),tau2tf(nt),tauntf(nt),gtf(nt), ...
                      thetatf(nt),phitf(nt),xtf(nt),a1,a2,rho,psi]);
    n1tf     = [n1t(2:nt);tem];
  end;
  n2t      = ( a2/a1* (1-phit)./(1-thetat).*((k2t./zt).^(phit*rho)./ ...
               ((k1t./zt).^(thetat*rho))).*n1t.^(1-rho*(1-thetat)) ).^ ...
               (ones(nt,1)./(1-rho*(1-phit)));
  n2tf     = ( a2/a1* (1-phitf)./(1-thetatf).*((k2tf./ztf).^(phitf*rho)./ ...
               ((k1tf./ztf).^(thetatf*rho))).*n1tf.^(1-rho*(1-thetatf))).^ ...
               (ones(nt,1)./(1-rho*(1-phitf)));
  y1t      = k1t.^thetat.*(zt.*n1t).^(1-thetat);
  y1tf     = k1tf.^thetatf.*(ztf.*n1tf).^(1-thetatf);
  y2t      = k2t.^phit  .*(zt.*n2t).^(1-phit);
  y2tf     = k2tf.^phitf.*(ztf.*n2tf).^(1-phitf);
  yt       = (a1*y1t.^rho+a2*y2t.^rho).^(1/rho);
  ytf      = (a1*y1tf.^rho+a2*y2tf.^rho).^(1/rho);
  ct       = yt-xt-gt;
  ctf      = ytf-xtf-gtf;
  p1t      = a1 *(yt./y1t).^(1-rho);
  p1tf     = a1 *(ytf./y1tf).^(1-rho);
  p2t      = a2 *(yt./y2t).^(1-rho);
  p2tf     = a2 *(ytf./y2tf).^(1-rho);
  r1t      = p1t.*thetat.*y1t./k1t;
  r1tf     = p1tf.*thetatf.*y1tf./k1tf;
  r2t      = p2t.*phit.*y2t./k2t;
  r2tf     = p2tf.*phitf.*y2tf./k2tf;
  wt       = p1t.*(1-thetat).*y1t./n1t; 
  wtf      = p1tf.*(1-thetatf).*y1tf./n1tf; 
  e1t      = (r1t-del1).*k1t;
  e1tf     = (r1tf-del1).*k1tf;
  ate1t    = (1-tau1t).*e1t;
  ate1tf   = (1-tau1tf).*e1tf;
  dt       = ate1t-k1tf+k1t;
  dtf      = ate1tf-k1tf2+k1tf;
  Retx     = (dtf+k1tf2)./k1tf-1;
  tem      = cumprod(1+Retx);
  mRetx(i) = (tem(length(tem))^(1/nyr)-1)*100;
  gnpt     = yt+tauct.*ct;
  dgnpt    = log(gnpt)-hptrend(log(gnpt),1600);
  dct      = log((1+tauct).*ct)-hptrend(log((1+tauct).*ct),1600);
  dxt      = log(xt)-hptrend(log(xt),1600);
  dgt      = log(gt)-hptrend(log(gt),1600);
  dnt      = log(n1t+n2t)-hptrend(log(n1t+n2t),1600);
  dx1t     = zeros(4,nyr); dx2t=dx1t;
  dx1t(:)  = x1t;  dx1t=sum(dx1t)';
  dx2t(:)  = x2t;  dx2t=sum(dx2t)';
  dx1t     = log(dx1t)-hptrend(log(dx1t),100);
  dx2t     = log(dx2t)-hptrend(log(dx2t),100);
  det      = log(ate1t)-hptrend(log(ate1t),1600);
  S(i,:)   = [std(dgnpt),std(dct),std(dxt),std(dgt),std(dnt),std(thetat)]*100;
  tem1     = corrcoef([thetat(1:nt-1),thetat(2:nt)]);
  tem2     = corrcoef([thetat,dgnpt]);
  tem3     = corrcoef([det,dct]);
  tem4     = corrcoef([dx1t,dx2t]);
  C(i,:)   = [tem1(1,2),tem2(1,2),tem3(1,2),tem4(1,2)];

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

disp(' ')
disp(' ')
disp(' ')
disp('MEANS OF RATIOS WITH RESPECT TO GNP')
gnp   = y+mtauc*c;
tab1  = [(1+mtauc)*c,x,mg,dep,gnp-dep, ...
          mtauc*c,y-dep,w*(n1+n2)+e1+e2,w*(n1+n2),e2,e1,mtau1*e1,ate1, ...
          ate1-re,re,p1*y1,dep1,w*n1,e1, ...
          (k1f+k2f)/4,k1f/4,k2f/4,x,x1,x2,dep,dep1,dep2, ...
          rev,mtauc*c,mtaun*w*(n1+n2),mtau1*e1,mtau2*e2,k1f/4]/gnp;
tab1  = [tab1,ate1*400/k1f];
str1  = ['  Components of GNP                                       \n';
         '    Private Consumption, (1+tau_c)*c                 %6.3f\n';
         '    Investment, x_1+x_2                              %6.3f\n';
         '    Government Purchases, g                          %6.3f\n';
         '                                                          \n';
         '  GNP --> Income                                          \n';
         '    - Depreciation, del_1 *k_1 + del_2 *k_2          %6.3f\n';
         '    = Net Product                                    %6.3f\n';
         '    - Consumption Taxes, tau_c*c                     %6.3f\n';
         '    = Income                                         %6.3f\n'; 
         '                                                          \n';
         '  Components of Income                               %6.3f\n';
         '    Wages, w*(n_1+n_2)                               %6.3f\n'; 
         '    Noncorporate Earnings, (r_2-del_2)*k_2           %6.3f\n'; 
         '    Corporate Earnings, (r_1-del_1)*k_1 = e_1        %6.3f\n'; 
         '      Tax, tau_1*e_1                                 %6.3f\n';
         '      Earnings after Tax, (1-tau_1)*e_1              %6.3f\n';
         '        Dividends, (1-tau_1)*e_1 -k_1''+k_1           %6.3f\n';
         '        Retained earnings, k_1''-k_1                  %6.3f\n'; 
         '                                                          \n';
         '  Corporate Output, p_1*y_1                          %6.3f\n';
         '    Depreciation, del_1*k_1                          %6.3f\n';
         '    Wages, w*n_1                                     %6.3f\n';
         '    Corporate Earnings, (r_1-del_1)*k_1              %6.3f\n';         
         '                                                          \n';
         '  Capital, end of year, (k_1''+k_2'')/4                %6.3f\n';
         '    Corporate, k_1''/4                                %6.3f\n';
         '    Noncorporate, k_2''/4                             %6.3f\n';
         '                                                          \n';
         '  Investment, x_1+x_2                                %6.3f\n';
         '    Corporate, x_1                                   %6.3f\n';
         '    Noncorporate, x_2                                %6.3f\n';  
         '                                                          \n';
         '  Depreciation, del_1*k_1 + del_2*k_2                %6.3f\n';
         '    Corporate, del_1*k_1                             %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 taxes, tau_1*e_1                       %6.3f\n';
         '    Noncorporate taxes, tau_2*e_2                    %6.3f\n';
         '                                                          \n';
         '  Value of Corporate Equities                        %6.3f\n';
         '    ==> After-tax earnings/Corporate value *400      %6.3f\n'];

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

disp('STANDARD DEVIATIONS OF LOGGED AND HP-FILTERED SERIES (%)')
str2   =['  GNP, (1+tau_c)*c+x+g                               %6.3f\n';
         '  Consumption, (1+tau_c)*c                           %6.3f\n';
         '  Investment, x                                      %6.3f\n';
         '  Government Purchases, g                            %6.3f\n';
         '  Hours of work, n                                   %6.3f\n'];
tab2   = [mean(S(:,1)),mean(S(:,2)),mean(S(:,3)),mean(S(:,4)),mean(S(:,5))];
disp(sprintf(str2',tab2))
disp(' ')

disp('CORRELATIONS FOR LOGGED AND HP-FILTERED VARIABLES')
str3   =['  Consumption & After-Tax earnings                   %6.3f\n';
         '  Corporate & Noncorporate Investment (annual)       %6.3f\n']; 
tab3   = [mean(C(:,3)),mean(C(:,4))];
disp(sprintf(str3',tab3))
disp(' ')

disp('STATISTICS FOR 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  = [mthet,mean(S(:,6)),mean(C(:,1)),mean(C(:,2))];
disp(sprintf(str4',tab4))
disp(' ')

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

