%SIM_MODEL    Use Vaughan's method to solve the neoclassical growth
%             model with input-output linkages:
%
%                max sum_t beta^t u(C[t]/N[t],L[t]/N[t])*N[t]
%
%             subj. to  sum_j Pj[t]Cj[t] <= sum_j Wj[t]Hj[t] 
%                         + sum_j (Dj[t]+Vj[t])Sj[t]-Vj[t]Sj[t+1]
%                         -Taxes+Transfers
%
%                    Taxes = tauc[t]*sum_j Pj[t]Cj[t] 
%                          + tauh[t]*sum_j Wj[t]Hj[t]
%                          + taud[t]*sum_j Dj[t]Sj[t]
%                                 - taux[t] X[t]-tauk K[t]}
%
%                   Firms: Yj[t] = F(Kj[t],Mj[t],Z[t]H[t])
%                   Resources: Yj[t] = Cj[t]+Xj[t]+sum_l Mjl[t]+Gj[t] 
%                   Utility: u(c,l) = {[c l^psi]^(1-alpha)-1}/(1-alpha)
%
%            with Z[t] = z[t](1+gz)^t and N[t]=(1+gn)^t.
%
%       

%        Ellen McGrattan, 12/15/13
%        Revised, ERM, 6/19/15

%---------------------------------------------------------------------
% PARAMETERS
%    

load ../MLEparam/fixedpar.dat
load ../MLEest/mleest.out
Theta0        = mleest;
T             = 10000;
nsim          = 1;
T             = 500;
nsim          = 1;

n             = fixedpar(1);
o             = ones(n,1);
alpha         = fixedpar(2);
beta          = fixedpar(3);
cbar          = fixedpar(4);
deli          = fixedpar(5:4+n);
delt          = fixedpar(5+n:4+2*n);
g             = fixedpar(5+2*n:4+3*n);
gamma         = zeros(n);
gamma(:)      = fixedpar(5+3*n:4+n*n+3*n);
gn            = fixedpar(5+n*n+3*n);
gz            = fixedpar(6+n*n+3*n);
nu            = zeros(n);
nu(:)         = fixedpar(7+n*n+3*n:6+2*n*n+3*n);
omega         = fixedpar(7+2*n*n+3*n:6+2*n*n+4*n);
p1            = fixedpar(7+2*n*n+4*n);
phi           = fixedpar(8+2*n*n+4*n:7+2*n*n+5*n);
psi           = fixedpar(8+2*n*n+5*n);
rho           = fixedpar(9+2*n*n+5*n);
tauc          = fixedpar(10+2*n*n+5*n);
taud          = fixedpar(11+2*n*n+5*n);
tauh          = fixedpar(12+2*n*n+5*n);
tauk          = fixedpar(13+2*n*n+5*n);
taup          = fixedpar(14+2*n*n+5*n);
taux          = fixedpar(15+2*n*n+5*n);
theta         = fixedpar(16+2*n*n+5*n:15+2*n*n+6*n);
z1            = fixedpar(16+2*n*n+6*n:15+2*n*n+7*n);
z2            = fixedpar(16+2*n*n+7*n:15+2*n*n+8*n);
zeta          = zeros(n);
zeta(:)       = fixedpar(16+2*n*n+8*n:15+3*n*n+8*n);
z             = fixedpar(16+3*n*n+8*n);
xguess        = fixedpar(17+3*n*n+8*n:16+3*n*n+11*n);
pvec          = fixedpar(17+3*n*n+11*n:16+3*n*n+12*n);
qvec          = fixedpar(17+3*n*n+12*n:16+3*n*n+13*n);
yvec          = fixedpar(17+3*n*n+13*n:16+3*n*n+14*n);
xvec          = fixedpar(17+3*n*n+14*n:16+3*n*n+15*n);
agghss        = fixedpar(17+3*n*n+15*n);

param         = [n;alpha;beta;cbar;deli;delt;g;gamma(:);gn;gz;nu(:);omega;
                 p1;phi;psi;rho;tauc;taud;tauh;tauk;taup;taux;theta;z1;z2;
                 zeta(:);z];

ns            = 3*n+7;
P             = zeros(ns);
Q             = zeros(ns);
ii            = 1:2*n;
P(ii,ii)      = .9*eye(2*n);
ii            = 1:n;
Q(ii,ii)      = .01*eye(n);
ii            = n+1:2*n;
Q(ii,ii)      = .0001*eye(n);

P(1,1)        = .995;
P(2,2)        = .995;
P(4,4)        = .995;
P(5,5)        = .995;
P(6,6)        = .995;
P(8,8)        = .995;
P(20,20)      = .995;
P(26,26)      = .995;
P(ns,ns)      = .995;
Q(ns,ns)      = 0.0104;
%
% Turn off cyclical variation in taul
%
P(48,48)      = .95;
Q(48,48)      = 0.00001;
wgt           = ones(n,1);

P(3,3)        = Theta0(1);
P(7,7)        = Theta0(2);
P(9,9)        = Theta0(3);
P(10,10)      = Theta0(4);
P(11,11)      = Theta0(5);
P(12,12)      = Theta0(6);
P(13,13)      = Theta0(7);
P(14,14)      = Theta0(8);
P(24,24)      = Theta0(9);
Q(1,1)        = Theta0(10);
Q(2,2)        = Theta0(11);
Q(3,3)        = Theta0(12);
Q(4,4)        = Theta0(13);
Q(5,5)        = Theta0(14);
Q(6,6)        = Theta0(15);
Q(7,7)        = Theta0(16);
Q(8,8)        = Theta0(17);
Q(9,9)        = Theta0(18);
Q(10,10)      = Theta0(19);
Q(11,11)      = Theta0(20);
Q(12,12)      = Theta0(21);
Q(13,13)      = Theta0(22);
Q(14,14)      = Theta0(23);
Q(20,20)      = Theta0(24);
Q(24,24)      = Theta0(25);
Q(26,26)      = Theta0(26);
Q(5,20)       = Theta0(27);
Q(9,24)       = Theta0(28);
Q(11,26)      = Theta0(29);
wgt           = ones(n,1);
wgt(1)        = Theta0(30);
wgt(2)        = Theta0(31);
wgt(3)        = Theta0(32);
wgt(4)        = Theta0(33);
wgt(5)        = Theta0(34);
wgt(6)        = Theta0(35);
wgt(7)        = Theta0(36);
wgt(8)        = Theta0(37);
wgt(9)        = Theta0(38);
wgt(10)       = Theta0(39);
wgt(11)       = Theta0(40);
wgt(12)       = Theta0(41);
wgt(13)       = Theta0(42);
wgt(14)       = Theta0(43);

Qhld=Q;
Q=.00001*Q;
Q(II,JJ)=Qhld(II,JJ);

%---------------------------------------------------------------------
% COMPUTE EQUILIBRIUM 
%    
%     X_{t+1} = A X_t + B S_t
%     Z_t     = C X_t + D S_t
%
%     where X_t = [ln kt_t, ln ki_t] 
%           Z_t = X_t
%           S_t = [ln z1_t,ln z2_t, ln g_t, tauc_t,taud_t,tauh_t, ...
%                  tauk_t,taup_t,taux_t, ln z_t]
%---------------------------------------------------------------------
%
% Step 1. Steady state 
%

x1         = ss_fullmodel(xguess,param);
o          = ones(n,1);
betg       = beta*(1+gz)^(-alpha);
gam        = sum(gamma)';
lshare     = 1-theta-phi-gam;
thet1      = 1-theta;
aggcnrm    = x1(1);
pvec       = [p1;x1(2:n)];
h2h1       = x1(n+1:2*n);
h          = x1(2*n+1:3*n);
pmvec      = gam.*prod(((pvec*o')./(gamma+1e-15)).^(gamma./(o*gam')))';
ptvec      = prod(((pvec*o')./(zeta+1e-15)).^zeta)';
aggp       = sum(omega.^rho.*pvec.^(1-rho))^(1/(1-rho));
cnrm       = (aggp*omega./pvec).^rho * aggcnrm;
c          = cnrm*cbar;
aggc       = (sum(omega.*c.^((rho-1)/rho)))^(rho/(rho-1));
aggl       = 1-sum(h);
m1y        = gam.*pvec./pmvec;
h1         = h./(1+h2h1);
h2         = h-h1;
qxipy      = h2h1;
xiy        = qxipy.^(1-phi).*(z2./z1).^lshare;
qp         = qxipy./xiy;
qvec       = qp.*pvec;
pivec      = prod(((qvec*o')./(nu+1e-15)).^nu)';
kiy        = betg*phi.*(qvec.*xiy+pvec)./((1-betg*(1-deli)).*pivec);
ret        = ((1+taux)*ptvec./pvec-betg*(1-(1-taup)*(delt+tauk) + ...
            (1-delt).*((1+taux)*ptvec./pvec-1)))/(betg*(1-taup));
k1ty       = theta./ret;
y          = (k1ty.^theta.*kiy.^phi.*m1y.^gam).^(o./lshare).*z1.*h1*z;
k1t        = k1ty.*y;
k2t        = h2h1.*k1t;
kt         = k1t+k2t;
m1         = m1y.*y;
xi         = xiy.*y;
ki         = kiy.*y;
qxi        = qvec.*xi;
qxip       = qxi./pvec;
prodxt     = ((1+gz)*(1+gn)-1+delt).*kt;
prodxi     = ((1+gz)*(1+gn)-1+deli).*ki;
m1mat      = gamma.*((o./pvec)*pvec').*(o*y');
m2mat      = gamma.*((o./pvec)*pvec').*(o*qxip');
xtmat      = zeta.*((o./pvec)*ptvec').*(o*prodxt');
ximat      = nu.*((o./qvec)*pivec').*(o*prodxi');
wage       = lshare.*pvec.*y./h1;
aggl       = 1-sum(h);
m2         = (prod(m2mat.^gamma)').^(o./gam);
ngdpnew    = sum(pvec.*y+qxi - pmvec.*(m1+m2));
rgdpnew    = ngdpnew;
by         = 1; 
idivcfti   = [ones(10,1);.69;ones(3,1);.85];
ngdpcur    = ngdpnew-qvec(11)*ximat(11,:)*idivcfti;
rgdpcur    = ngdpcur;
ngdpold    = ngdpnew-sum(sum((qvec*o').*ximat));
rgdpold    = ngdpold;
tfpold     = rgdpold/(sum(pvec.*kt)^.33*sum(h)^.67);
tfpcur     = rgdpcur/(sum(pvec.*kt)^.33*sum(h)^.67);
tfpnew     = rgdpnew/(sum(pvec.*kt)^.33*sum(h)^.67);
lwold      = (psi*(1+tauc)*aggp*aggc/aggl)/(.67*(1-tauh)*ngdpold/(1-aggl));
lwcur      = (psi*(1+tauc)*aggp*aggc/aggl)/(.67*(1-tauh)*ngdpcur/(1-aggl));
lwnew      = (psi*(1+tauc)*aggp*aggc/aggl)/(.67*(1-tauh)*ngdpnew/(1-aggl));
yvec       = y;
xvec       = xi;
Sbar       = [log(z1);log(z2);log(g);tauc;taud;tauh;tauk;taup;taux;log(z)];

%
% Step 2. Linearize first-order conditions
%
m1matplus  = max(m1mat,1e-10);
m2matplus  = max(m2mat,1e-10);
v          = [log(kt);log(ki);log(kt);log(ki);log(kt);log(ki);Sbar;Sbar;
              log(c);log(h1);log(h2);log(m1matplus(:));log(m2matplus(:));
              log(c);log(h1);log(h2);log(m1matplus(:));log(m2matplus(:))];
del        = max(abs(v)*1e-4,1e-8);
nres       = 8*n+4*n*n;
nv         = length(v);
[res,xxs,xxxs] = ...
             res_fullmodel(v,[param;pvec;qvec;yvec;xvec;1-aggl;
             h2([5,9,11]);wgt]);
dR         = zeros(nres,nv);
dV         = zeros(length(xxs(:)),nv);
dVV        = zeros(length(xxxs(:)),nv);

for i=1:nv;
  vp           = v;
  vm           = v;
  vp(i)        = vp(i)+del(i);
  vm(i)        = vm(i)-del(i);
  [rp,xp,xxp]  = res_fullmodel(vp,[param;pvec;qvec;yvec;xvec;1-aggl;
                 h2([5,9,11]);wgt]);
  [rm,xm,xxm]  = res_fullmodel(vm,[param;pvec;qvec;yvec;xvec;1-aggl;
                 h2([5,9,11]);wgt]);
  dR(:,i)      = (rp-rm)/(2*del(i));
  dV(:,i)      = (xp-xm)/(2*del(i));
  dVV(:,i)     = (xxp-xxm)/(2*del(i));
end;
i1         = 1:6*n+4*n*n;
i2         = 6*n+4*n*n+1:8*n+4*n*n;
j1         = 1:6*n+2*ns;
j2         = 6*n+2*ns+1:nv;
Coef       = -inv(dR(i1,j2))*dR(i1,j1);
dRnew      = dR(i2,j1)+dR(i2,j2)* Coef;
dVnew      = dV(:,j1)+dV(:,j2)* Coef;
dVVnew     = dVV(:,j1)+dVV(:,j2)* Coef;

%
% Step 2. Find deterministic part of solution
%
%         Residuals are A1*[X_{t+1}] + A2*[X_t]  + shock terms = 0
%                          [Z_{t+1}]      [Z_t]
%
%                           X_t  = [kt[t],ki[t]], 
%                           Z_t  = [X_t,h1[t],h2[t],y[t]]
%

nX         = 2*n;
nZ         = 2*n;
nXZ        = 4*n;

a1         = dRnew(:,1:2*n);
a2         = dRnew(:,2*n+1:4*n);
a3         = dRnew(:,4*n+1:6*n);
a4         = dRnew(:,6*n+1:6*n+ns);
a5         = dRnew(:,6*n+ns+1:6*n+2*ns);

A1         = [eye(2*n),zeros(2*n);
              zeros(2*n),a3];
A2         = [zeros(2*n),-eye(2*n);
              a1,a2];

[V,L]      = eig(A2,-A1);
L          = diag(L);
[D,I]      = sort(abs(L));
if (~( (D(nX) <= 1) & (D(nX+1) > 1) ) )
  i      = sum(D<=1);
  disp(sprintf( ...
  'WARNING: program assumes %g eigenvalues inside unit circle, found %g',nX,i));
end;
V11        = V(1:nX,I(1:nX));
V21        = V(nX+1:nXZ,I(1:nX));
W          = real([V11*diag(L(I(1:nX)))/V11;
                   V21/V11]);
A          = W(1:2*n,:);
C          = W(2*n+1:4*n,:);

%
% Step 3. Find the stochastic part of the solution (ie., matrices B,D)
%
%         Residuals (k in logs): 
%
%           E_t{ a1 k[t]+ a2 k[t+1] + a3 k[t+2] + 
%                         a4 S[t] + a5 S[t+1] } = 0,
%
%           where k =[ln(kt);ln(ki)]
%                 S =[ln(z1),ln(z2),ln(g),tauc,taud,tauh,tauk,taup,taux,ln(z)] 
%                     (3n+7 x 1)
%
%           ==>  (a1+a2*A+a3*A*A)*k
%              + (a2*B+a3*A*B+a3*B*P+a4+a5*P)*S = 0
%
%           Note: need vec(A+B)=vec(A)+vec(B), vec(ABC)=(C' kron A)*vec(B)
%
%

resk       = a1+a2*A+a3*A*A;

if (any(abs(resk)>1e-8));
  disp('Check deterministic part of the solution')
end;

I          = eye(ns);
tem        =-[kron(I,a2+a3*A)+kron(P',a3)]\vec(a4+a5*P);
B          = zeros(2*n,ns);
B(:)       = tem(1:2*n*ns);
D          = B;

ii         = [1:3*n+2*n*n];   % indices for c[t],h1[t],h2[t],m1mat[t],m2mat[t]  
C2         = Coef(ii,1:2*n)+Coef(ii,2*n+1:4*n)*A;
D2         = Coef(ii,6*n+1:6*n+ns)+Coef(ii,2*n+1:4*n)*B;

%
% Add additional variables to vector of observables:
%
C2         = [C2;dVnew(:,1:2*n)+dVnew(:,2*n+1:4*n)*A;
                 dVVnew(:,1:2*n)+dVVnew(:,2*n+1:4*n)*A]; 
D2         = [D2;dVnew(:,6*n+1:6*n+ns)+dVnew(:,2*n+1:4*n)*B;
                 dVVnew(:,6*n+1:6*n+ns)+dVVnew(:,2*n+1:4*n)*B];

%
% Step 4. Simulate time series for x=[kt,ki,Sbar], u=eps,
%                                  y=[c,h1,h2,m1mat,m2mat], all in logs
%
%           x[n+1] = [A,  B,  A0] x[n] + [ 0 ] u[n]
%                    [0   P,  P0]        [ Q ] 
%                    [0   0,   1]        [ 0 ]
%  
%                        (G)              (H)
%
%           y[n]   = [C,  D,  C0] x[n] 
%
%                        (J)
%
A0         = (eye(2*n)-A)*log([kt;ki])-B*Sbar;
C0         = [log([c;h1;h2;m1matplus(:);m2matplus(:)]);xxs;xxxs]- ...
              C2*log([kt;ki])-D2*Sbar;
nc         = length(C0);
P0         = (I-P)*Sbar;
G          = [A,B,A0;zeros(3*n+7,2*n),P,P0;zeros(1,5*n+7),1];
H          = [zeros(2*n,3*n+7);Q;zeros(1,3*n+7)];
J          = [C2,D2,C0];
J0         = J(1:495,:);
J1         = J(496:513,:);
J2         = J(514:628,:);
X0         = [log(kt);log(ki);Sbar;1];
m1matt     = zeros(T,n*n);
m2matt     = zeros(T,n*n);
%rng(3,'twister')

%rng(1,'twister')

hp         = zeros(nsim,17);
lw         = zeros(nsim,4);
ot         = ones(T,1);
tem        = zeros(n);
for s=1:nsim
  eps      = randn(T+1,ns);
  [Z,X]    = dlsim(G,H,J,zeros(nc,3*n+7),eps,X0); 
  ct       = exp(Z(1:T,1:n));
  h1t      = exp(Z(1:T,n+1:2*n));
  h2t      = exp(Z(1:T,2*n+1:3*n));
  m1matt(:)= exp(Z(1:T,3*n+1:3*n+n*n));
  m2matt(:)= exp(Z(1:T,3*n+n*n+1:3*n+2*n*n));
  ktt      = exp(X(1:T,1:n));
  kttp     = exp(X(2:T+1,1:n));
  kit      = exp(X(1:T,n+1:2*n));
  kitp     = exp(X(2:T+1,n+1:2*n));
  zvar     = Z(1:T,3*n+2*n*n+1:nc);
  z1t      = exp(X(1:T,2*n+1:3*n));
  z2t      = exp(X(1:T,3*n+1:4*n));
  gt       = exp(X(1:T,4*n+1:5*n));
  tauct    = X(1:T,5*n+1);
  taudt    = X(1:T,5*n+2);
  tauht    = X(1:T,5*n+3);
  taukt    = X(1:T,5*n+4);
  taupt    = X(1:T,5*n+5);
  tauxt    = X(1:T,5*n+6);
  zt       = exp(X(1:T,5*n+7));
  ht       = h1t+h2t;
  qxipyt   = h2t./h1t;
  k1tt     = ktt./(ot*o'+qxipyt);
  k2tt     = ktt-k1tt;
  pt       = p1*(ct./(ct(:,1)*o')).^(-1/rho) .*(ot*omega')/omega(1);
  for t=1:T;
    tem(:)   = m1matt(t,:)';
    m1t(t,:) = ((prod(tem.^gamma)').^(o./gam))';
    tem(:)   = m2matt(t,:)';
    m2t(t,:) = ((prod(tem.^gamma)').^(o./gam))';
  end;
  zzt      = (zt*o').^(ot*wgt');
  zzt(:,n) = ot;
  yt       = k1tt.^(ot*theta').*kit.^(ot*phi').*m1t.^(ot*gam').* ...
             (zzt.*z1t.*h1t).^(ot*lshare');
  xit      = k2tt.^(ot*theta').*kit.^(ot*phi').*m2t.^(ot*gam').* ...
             (zzt.*z2t.*h2t).^(ot*lshare');
  prodxtt  = (1+gz)*(1+gn)*kttp-(ot*(1-delt)').*ktt;
  prodxit  = (1+gz)*(1+gn)*kitp-(ot*(1-deli)').*kit;
  qxipt    = qxipyt.*yt;
  qxit     = qxipt.*pt;
  qt       = qxit./xit;
  for t=1:T;
    pmt(t,:) = (gam.*prod(((pt(t,:)'*o')./gamma).^(gamma./(o*gam')))')';
    ptt(t,:) = prod(((pt(t,:)'*o')./(zeta+1e-15)).^zeta);
    pit(t,:) = prod(((qt(t,:)'*o')./(nu+1e-15)).^nu);
    tem      = zeta.*((o./pt(t,:)')*ptt(t,:)).*(o*prodxtt(t,:));
    xtt(t,:) = sum(tem');
    tem      = nu.*((o./qt(t,:)')*pit(t,:)).*(o*prodxit(t,:));
    xit2(t,:)= sum(tem');
    xi_row_11(t,:) = tem(11,:);
  end;
  if n>1;
    aggct  = (sum( ((ot*omega').*ct.^((rho-1)/rho))' )').^(rho/(rho-1));
    agglt  = 1-sum(ht')';
    aggpt  = (sum(((ot*omega').^rho.*pt.^(1-rho))')').^(1/(1-rho));
  else
    aggct  = (omega*ct.^((rho-1)/rho)).^(rho/(rho-1));
    agglt  = 1-ht;
    aggpt  = (omega^rho*pt.^(1-rho)).^(1/(1-rho));
  end;
  agght    = 1-agglt;
  % 
  % GDP and the labor wedge
  %
  if (n>1); 
    ngdpnewt = sum( (pt.*yt+qt.*xit - pmt.*m1t -pmt.*m2t)')';
    rgdpnewt = sum( ((ot*pt(by,:)).*yt+(ot*qt(by,:)).*xit - ...
               (ot*pmt(by,:)).*(m1t+m2t))')';
    rgdpnewt2= fisher([pt,qt,-pmt],[yt,xit,m1t+m2t],by,1);
    rbusrect = sum( ((ot*pt(by,:)).*yt+(ot*qt(by,:)).*xit)')';
    rbusrect = sum( ((ot*pvec').*yt+(ot*qvec').*xit)')';
    riinvt   = sum( ((ot*qvec').*xit)')';
    ngdpoldt = sum( (pt.*yt - pmt.*m1t -pmt.*m2t)')';
    rgdpoldt = sum( ((ot*pt(by,:)).*yt - (ot*pmt(by,:)).*(m1t+m2t))')';
    naggpyt  = sum( (pt.*yt)' )';
    raggpyt  = sum( ((ot*pt(by,:)).*yt)' )';
    naggpxtt = sum( (pt.*xtt)' )';
    raggpxtt = sum( ((ot*pt(by,:)).*xtt)' )';
    naggpktt = sum( (pt.*ktt)' )';
    raggpktt = sum( ((ot*pt(by,:)).*ktt)' )';
    naggqxit = sum( (qxit)' )';
    raggqxit = sum( ((ot*qt(by,:)).*xit)' )';
    naggqkit = sum( (qt.*kit)' )';
    raggqkit = sum( ((ot*qt(by,:)).*kit)' )';
    naggpgt  = sum( ((ot*pt(by,:)).*gt)' )';
    raggpgt  = sum( (pt.*gt)' )';
    naggpct  = sum( (pt.*ct)' )'; 
    raggpct  = sum( ((ot*pt(by,:)).*ct)' )'; 
    aggh1t   = sum(h1t')';
    aggh2t   = sum(h2t')';
  else
    ngdpnewt = pt.*yt +qt.*xit - pmt.*m1t -pmt.*m2t;
    rgdpnewt = (ot*pt(by,:)).*yt+(ot*qt(by,:)).*xit-(ot.*pmt(by,:)).*(m1t+m2t);
    ngdpoldt = pt.*yt - pmt.*m1t -pmt.*m2t;
    rgdpoldt = (ot*pt(by,:)).*yt-(ot.*pmt(by,:)).*(m1t+m2t);
    naggpyt  = pt.*yt;
    raggpyt  = (ot*pt(by,:)).*yt;
    naggpxtt = pt.*xtt;
    raggpxtt = (ot*pt(by,:)).*xtt;
    naggpktt = pt.*ktt;
    raggpktt = (ot*pt(by,:)).*ktt;
    naggqxit = qxit;
    raggqxit = (ot*qt(by,:)).*xit;
    naggqkit = qt.*kit;
    raggqkit = (ot*qt(by,:)).*kit;
    naggpgt  = pt.*gt;
    raggpgt  = (ot*pt(by,:)).*gt;
    naggpct  = pt.*ct;
    raggpct  = (ot*pt(by,:)).*ct;
    aggh1t   = h1t;
    aggh2t   = h2t;
  end;
  rgdpochkt = fisher([pt,-pmt],[yt,m1t+m2t],by,1);
  rgot      = fisher([pt,qt],[yt,xit],by,1);
  ngdpcurt  = ngdpnewt-qt(:,11).* (xi_row_11*idivcfti); 
  rgdpcurt  = rgdpnewt-qt(by,11)* (xi_row_11*idivcfti); 
  rtfpoldt  = rgdpoldt./(raggpktt.^.33.*agght.^.67);
  rtfpcurt  = rgdpcurt./(raggpktt.^.33.*agght.^.67);
  rtfpnewt  = rgdpnewt./(raggpktt.^.33.*agght.^.67);
  ntfpoldt  = ngdpoldt./(naggpktt.^.33.*agght.^.67);
  ntfpcurt  = ngdpcurt./(naggpktt.^.33.*agght.^.67);
  ntfpnewt  = ngdpnewt./(naggpktt.^.33.*agght.^.67);
  lwoldt    = (psi*(1+tauct).*naggpct./agglt)./(.67*(1-tauht).* ...
              ngdpoldt./(1-agglt));
  lwcurt    = (psi*(1+tauct).*naggpct./agglt)./(.67*(1-tauht).* ...
              ngdpcurt./(1-agglt));
  lwnewt    = (psi*(1+tauct).*naggpct./agglt)./(.67*(1-tauht).* ...
              ngdpnewt./(1-agglt));
  defoldt   = ngdpoldt./rgdpoldt;

  for i=1:n;
  lwt(:,i)  = (psi*(1+tauct).*naggpct./agglt)./(lshare(i)*(1-tauht).* ...
               pt(:,i).*yt(:,i)./h1t(:,i));
  end;

  lyt       = log(raggpyt)-hptrend(log(raggpyt),1600);
  lct       = log(raggpct)-hptrend(log(raggpct),1600);
  lxtt      = log(raggpxtt)-hptrend(log(raggpxtt),1600);
  lxit      = log(raggqxit)-hptrend(log(raggqxit),1600); 
  lht       = log(agght)-hptrend(log(agght),1600);
  lh1t      = log(aggh1t)-hptrend(log(aggh1t),1600);
  lh2t      = log(aggh2t)-hptrend(log(aggh2t),1600);
  llwt      = log(lwoldt)-hptrend(log(lwoldt),1600);
  llwnt     = log(lwnewt)-hptrend(log(lwnewt),1600);
  clyct     = corrcoef([lyt,lct]);
  clyxtt    = corrcoef([lyt,lxtt]);
  clyxit    = corrcoef([lyt,lxit]);
  clyht     = corrcoef([lyt,lht]);
  clyh1t    = corrcoef([lyt,lh1t]);
  clyh2t    = corrcoef([lyt,lh2t]);
  clhlwt    = corrcoef([lht,llwt]);
  clhlwnt   = corrcoef([lht,llwnt]);
  hp(s,:)   = [[std(lyt),std(lct),std(lxtt),std(lxit),std(lht), ...
                std(lh1t),std(lh2t),std(llwt),std(llwnt)]*100, ...
               clyct(1,2),clyxtt(1,2),clyxit(1,2),clyht(1,2), ...
               clyh1t(1,2),clyh2t(1,2),clhlwt(1,2),clhlwnt(1,2)];

  lwvars    = log([lwoldt,agght,raggqxit,raggpxtt,lwnewt]);
  lwtem     = corrcoef(lwvars);
  lw(s,:)   = lwtem(1,2:5);

end;
if nsim>1;
  hptab     = [mean(hp)',std(hp)'];
  lwtab     = [mean(lw)',std(lw)'];
else
  hptab     = [hp',0*hp'];
  lwtab     = [lw',0*lw'];
end;

%
% Step 5. Print results
%
clc
disp('MULTI-SECTOR MODEL WITH INTANGIBLES, ELASTIC LABOR & DISTORTIONS ')
disp('---------------------------------------------------------------- ')
disp(' ')
disp('  Logged and HP-filtered statistics:')
disp(sprintf('    std(gross output)               %7.3f (%7.3f)',hptab(1,:)))
disp(sprintf('    std(consumption)                %7.3f (%7.3f)',hptab(2,:)))
disp(sprintf('    std(tangible inv.)              %7.3f (%7.3f)',hptab(3,:)))
disp(sprintf('    std(intangible inv.)            %7.3f (%7.3f)',hptab(4,:)))
disp(sprintf('    std(hours)                      %7.3f (%7.3f)',hptab(5,:)))
disp(sprintf('    std(hours, g&s)                 %7.3f (%7.3f)',hptab(6,:)))
disp(sprintf('    std(hours, intang)              %7.3f (%7.3f)',hptab(7,:)))
disp(sprintf('    std(labor wedge)                %7.3f (%7.3f)',hptab(8,:)))
disp(sprintf('    std(labor wedge, new)           %7.3f (%7.3f)',hptab(9,:)))
disp(sprintf('    cor(gross output,consumption)   %7.3f (%7.3f)',hptab(10,:)))
disp(sprintf('    cor(gross output,tangible inv)  %7.3f (%7.3f)',hptab(11,:)))
disp(sprintf('    cor(gross output,intang inv)    %7.3f (%7.3f)',hptab(12,:)))
disp(sprintf('    cor(gross output,hours)         %7.3f (%7.3f)',hptab(13,:)))
disp(sprintf('    cor(gross output,hours g&s)     %7.3f (%7.3f)',hptab(14,:)))
disp(sprintf('    cor(gross output,hours intang)  %7.3f (%7.3f)',hptab(15,:)))
disp(sprintf('    cor(labor wedge,hours)          %7.3f (%7.3f)',hptab(16,:)))
disp(sprintf('    cor(labor wedge new, hours)     %7.3f (%7.3f)',hptab(17,:)))
disp(' ')
disp('  Logged and unfiltered statistics:')
disp(sprintf('    cor(labor wedge,hours)          %7.3f (%7.3f)',lwtab(1,:)))
disp(sprintf('    cor(labor wedge,intangible inv) %7.3f (%7.3f)',lwtab(2,:)))
disp(sprintf('    cor(labor wedge,tangible inv)   %7.3f (%7.3f)',lwtab(3,:)))
disp(sprintf('    cor(labor wedge,labor wedge new)%7.3f (%7.3f)',lwtab(4,:)))
disp('---------------------------------------------------------------- ')
disp(sprintf(' Note: %4.0f simulations of length %4.0f',[nsim,T]))


return

disp(' Time series statistics:')
disp(' ')
disp(sprintf('    corr(hours,labor wedge, old)    %7.3f (%7.3f)',tab0(1,:)))
disp(sprintf('    corr(hours,labor wedge, new)    %7.3f (%7.3f)',tab0(2,:)))
disp(sprintf('    corr(hours,intangible inv)      %7.3f (%7.3f)',tab0(3,:)))
disp(' ')
disp(' Time series statistics (logged and detrended variables): ')
disp(' ')
disp(sprintf('    std(gdp current)        %7.3f (%7.3f)',tab(10,:)))
disp(sprintf('    std(gdp new)            %7.3f (%7.3f)',tab(11,:)))
disp(sprintf('    std(nom gdp old)        %7.3f (%7.3f)',tab(12,:)))
disp(sprintf('    std(nom gdp current)    %7.3f (%7.3f)',tab(13,:)))
disp(sprintf('    std(nom gdp new)        %7.3f (%7.3f)',tab(14,:)))
disp(sprintf('    std(labor wedge old)    %7.3f (%7.3f)',tab(15,:)))
disp(sprintf('    std(labor wedge cur)    %7.3f (%7.3f)',tab(16,:)))
disp(sprintf('    std(labor wedge new)    %7.3f (%7.3f)',tab(17,:)))
disp(sprintf('    std(tfp old)            %7.3f (%7.3f)',tab(18,:)))
disp(sprintf('    std(tfp cur)            %7.3f (%7.3f)',tab(19,:)))
disp(sprintf('    std(tfp new)            %7.3f (%7.3f)',tab(20,:)))
disp(sprintf('    cor(wedge, gdp old)     %7.3f (%7.3f)',tab(21,:)))
disp(sprintf('    cor(wedge, gdp cur)     %7.3f (%7.3f)',tab(22,:)))
disp(sprintf('    cor(wedge, gdp new)     %7.3f (%7.3f)',tab(23,:)))
disp(sprintf('    cor(hours, gdp old)     %7.3f (%7.3f)',tab(24,:)))
disp(sprintf('    cor(hours, gdp cur)     %7.3f (%7.3f)',tab(25,:)))
disp(sprintf('    cor(hours, gdp new)     %7.3f (%7.3f)',tab(26,:)))
disp(sprintf('    cor(tfp, gdp old)       %7.3f (%7.3f)',tab(27,:)))
disp(sprintf('    cor(tfp, gdp cur)       %7.3f (%7.3f)',tab(28,:)))
disp(sprintf('    cor(tfp, gdp new)       %7.3f (%7.3f)',tab(29,:)))



