This is the data file for Yi Wen's article "Granger Causality and Real Business Cycles" in the
May/June 2007 issue of the Federal Reserve Bank of St. Louis REVIEW.


0. US data files: 								.dat ; .txt

1. Rats files:
Granger test for Updated US data:						Granger_new#.prg
Granger test for artificial data generated from models:				Granger#.prg

2. Gauss files:
Standard RBC Model:								kpr1_rbc.gss
RBC model with sequential information structure:				Wen.gss
Scale economt model:								Granger.gss (use this model to generate figure 1)

3. outputs generated from Gauss files are stored in .out files.

*********************************************************************************************************************************
*********************************************************************************************************************************
*********************************************************************************************************************************

RATS Programs

*********************************************************************************************************************************
*	Granger_New1.prg													*
*********************************************************************************************************************************

cal 1947 1 4
all 2006:1
open data c:\rats\yi\granger\rev2006\gdpq1947.txt
data(for=free,org=obs) / gdpq
open data c:\rats\yi\granger\rev2006\gcq1947.txt
data(for=free,org=obs) / gcq
open data c:\rats\yi\granger\rev2006\gcnq1947.txt
data(for=free,org=obs) / gcnq
open data c:\rats\yi\granger\rev2006\gcsq1947.txt
data(for=free,org=obs) / gcsq
open data c:\rats\yi\granger\rev2006\ginq1947.txt
data(for=free,org=obs) / ginq
open data c:\rats\yi\granger\rev2006\gnxq1947.txt
data(for=free,org=obs) / gnxq
open data c:\rats\yi\granger\rev2006\ggq1947.txt
data(for=free,org=obs) / ggq
open data c:\rats\yi\granger\rev2006\givq1947.txt
data(for=free,org=obs) / givq

set y = gdpq
set c = gcnq+gcsq
set i1 = ginq
set iv = givq

set trend = 1
set tt = trend**2
linreg(noprint) y / resy
#constant trend tt
set growth = y - resy
set iv = iv/growth

set y = gdpq-givq

set dy = log(y(t)) - log(y(t-1))
set dc = log(c(t)) - log(c(t-1))
set di = log(i1(t)) - log(i1(t-1))
set div = iv(t) - iv(t-1)

smpl 1947:1 2006:1


graph 3
#dy
#dc
#di


linreg dy
#constant trend dy{1 to 2}

linreg dy
#constant trend dy{1 to 2} dc{1}

linreg dy
#constant trend dy{1 to 2} di{1}

linreg dy
#constant trend dy{1 to 2} div{1}

linreg dc
#constant trend dc{1 to 2}

linreg dc
#constant trend dc{1 to 2} dy{1}

linreg dc
#constant trend dc{1 to 2} di{1}

linreg dc
#constant trend dc{1 to 2} div{1}

linreg di
#constant trend di{1 to 2}

linreg di
#constant trend di{1 to 2} dc{1}

linreg di
#constant trend di{1 to 2} dy{1}

linreg di
#constant trend di{1 to 2} div{1}

linreg div
#constant trend div{1 to 2}

linreg div
#constant trend div{1 to 2} dy{1}

linreg div
#constant trend div{1 to 2} dc{1}

linreg div
#constant trend div{1 to 2} di{1}


*********************************************************************************************************************************
*	Granger_New2.prg													*
*********************************************************************************************************************************

cal 1947 1 4
all 2006:1
open data c:\rats\yi\granger\rev2006\gdpq1947.txt
data(for=free,org=obs) / gdpq
open data c:\rats\yi\granger\rev2006\gcq1947.txt
data(for=free,org=obs) / gcq
open data c:\rats\yi\granger\rev2006\gcnq1947.txt
data(for=free,org=obs) / gcnq
open data c:\rats\yi\granger\rev2006\gcsq1947.txt
data(for=free,org=obs) / gcsq
open data c:\rats\yi\granger\rev2006\ginq1947.txt
data(for=free,org=obs) / ginq
open data c:\rats\yi\granger\rev2006\gnxq1947.txt
data(for=free,org=obs) / gnxq
open data c:\rats\yi\granger\rev2006\ggq1947.txt
data(for=free,org=obs) / ggq
open data c:\rats\yi\granger\rev2006\givq1947.txt
data(for=free,org=obs) / givq

set y = gdpq
set c = gcnq+gcsq
set i1 = ginq
set iv = givq

set trend = 1
set tt = trend**2
linreg(noprint) y / resy
#constant trend tt
set growth = y - resy
set iv = iv/growth

set y = gdpq-givq

set dy = log(y(t)) - log(y(t-1))
set dc = log(c(t)) - log(c(t-1))
set di = log(i1(t)) - log(i1(t-1))
set div = iv(t) - iv(t-1)

smpl 1947:1 2006:1


graph 3
#dy
#dc
#di


linreg dy
#constant trend dy{1 to 2}

linreg dy
#constant trend dy{1 to 2} dc{1}

linreg dy
#constant trend dy{1 to 2} di{1}

linreg dy
#constant trend dy{1 to 2} div{1}

linreg dc
#constant trend dc{1 to 2}

linreg dc
#constant trend dc{1 to 2} dy{1}

linreg dc
#constant trend dc{1 to 2} di{1}

linreg dc
#constant trend dc{1 to 2} div{1}

linreg di
#constant trend di{1 to 2}

linreg di
#constant trend di{1 to 2} dc{1}

linreg di
#constant trend di{1 to 2} dy{1}

linreg di
#constant trend di{1 to 2} div{1}

linreg div
#constant trend div{1 to 2}

linreg div
#constant trend div{1 to 2} dy{1}

linreg div
#constant trend div{1 to 2} dc{1}

linreg div
#constant trend div{1 to 2} di{1}


*********************************************************************************************************************************
*	Granger1.prg														*
*********************************************************************************************************************************

cal 1947 1 4
all 2006:1
open data c:\rats\yi\granger\rev2006\gdpq1947.txt
data(for=free,org=obs) / gdpq
open data c:\rats\yi\granger\rev2006\gcq1947.txt
data(for=free,org=obs) / gcq
open data c:\rats\yi\granger\rev2006\gcnq1947.txt
data(for=free,org=obs) / gcnq
open data c:\rats\yi\granger\rev2006\gcsq1947.txt
data(for=free,org=obs) / gcsq
open data c:\rats\yi\granger\rev2006\ginq1947.txt
data(for=free,org=obs) / ginq
open data c:\rats\yi\granger\rev2006\gnxq1947.txt
data(for=free,org=obs) / gnxq
open data c:\rats\yi\granger\rev2006\ggq1947.txt
data(for=free,org=obs) / ggq
open data c:\rats\yi\granger\rev2006\givq1947.txt
data(for=free,org=obs) / givq

set y = gdpq
set c = gcnq+gcsq
set i1 = ginq
set iv = givq

set trend = 1
set tt = trend**2
linreg(noprint) y / resy
#constant trend tt
set growth = y - resy
set iv = iv/growth

set y = gdpq-givq

set dy = log(y(t)) - log(y(t-1))
set dc = log(c(t)) - log(c(t-1))
set di = log(i1(t)) - log(i1(t-1))
set div = iv(t) - iv(t-1)

smpl 1947:1 2006:1


graph 3
#dy
#dc
#di


linreg dy
#constant trend dy{1 to 2}

linreg dy
#constant trend dy{1 to 2} dc{1}

linreg dy
#constant trend dy{1 to 2} di{1}

linreg dy
#constant trend dy{1 to 2} div{1}

linreg dc
#constant trend dc{1 to 2}

linreg dc
#constant trend dc{1 to 2} dy{1}

linreg dc
#constant trend dc{1 to 2} di{1}

linreg dc
#constant trend dc{1 to 2} div{1}

linreg di
#constant trend di{1 to 2}

linreg di
#constant trend di{1 to 2} dc{1}

linreg di
#constant trend di{1 to 2} dy{1}

linreg di
#constant trend di{1 to 2} div{1}

linreg div
#constant trend div{1 to 2}

linreg div
#constant trend div{1 to 2} dy{1}

linreg div
#constant trend div{1 to 2} dc{1}

linreg div
#constant trend div{1 to 2} di{1}

*********************************************************************************************************************************
*	Granger2.prg														*
*********************************************************************************************************************************

cal 1947 1 4
all 2006:1
open data c:\rats\yi\granger\rev2006\gdpq1947.txt
data(for=free,org=obs) / gdpq
open data c:\rats\yi\granger\rev2006\gcq1947.txt
data(for=free,org=obs) / gcq
open data c:\rats\yi\granger\rev2006\gcnq1947.txt
data(for=free,org=obs) / gcnq
open data c:\rats\yi\granger\rev2006\gcsq1947.txt
data(for=free,org=obs) / gcsq
open data c:\rats\yi\granger\rev2006\ginq1947.txt
data(for=free,org=obs) / ginq
open data c:\rats\yi\granger\rev2006\gnxq1947.txt
data(for=free,org=obs) / gnxq
open data c:\rats\yi\granger\rev2006\ggq1947.txt
data(for=free,org=obs) / ggq
open data c:\rats\yi\granger\rev2006\givq1947.txt
data(for=free,org=obs) / givq

set y = gdpq
set c = gcnq+gcsq
set i1 = ginq
set iv = givq

set trend = 1
set tt = trend**2
linreg(noprint) y / resy
#constant trend tt
set growth = y - resy
set iv = iv/growth

set y = gdpq-givq

set dy = log(y(t)) - log(y(t-1))
set dc = log(c(t)) - log(c(t-1))
set di = log(i1(t)) - log(i1(t-1))
set div = iv(t) - iv(t-1)

smpl 1947:1 2006:1


graph 3
#dy
#dc
#di


linreg dy
#constant trend dy{1 to 2}

linreg dy
#constant trend dy{1 to 2} dc{1}

linreg dy
#constant trend dy{1 to 2} di{1}

linreg dy
#constant trend dy{1 to 2} div{1}

linreg dc
#constant trend dc{1 to 2}

linreg dc
#constant trend dc{1 to 2} dy{1}

linreg dc
#constant trend dc{1 to 2} di{1}

linreg dc
#constant trend dc{1 to 2} div{1}

linreg di
#constant trend di{1 to 2}

linreg di
#constant trend di{1 to 2} dc{1}

linreg di
#constant trend di{1 to 2} dy{1}

linreg di
#constant trend di{1 to 2} div{1}

linreg div
#constant trend div{1 to 2}

linreg div
#constant trend div{1 to 2} dy{1}

linreg div
#constant trend div{1 to 2} dc{1}

linreg div
#constant trend div{1 to 2} di{1}


*********************************************************************************************************************************
*********************************************************************************************************************************
*********************************************************************************************************************************

Gauss Programs

*********************************************************************************************************************************
*	Kpr1_rbc.gss														*
*********************************************************************************************************************************

/*
   kpr1_rbc.gss, YW, 5/16/95 -- Standard KPR Model

 -- AR(2) PERIODIC TECHNOLOGY SHOCKS:
    A(t) = r1 A(t-1) + r2 A(t-2) + e(t)
*/

@ ECONOMIC PARAMETER VALUES @

alp = 0.3;       @ Capital's share @
deltak=.025;      @ Depreciation Rate for capital @
gama=0.0000025;    @ > 0 singularity occurs at 0, Elasiticity of labor supply @
gammax=1.00;
r=.065/4;
beta=gammax/(1+r);
beta = 0.99;

g = 1;          @ std of tech shock @
rho1 = 1.6;     @ Serial Correlation of tech. shock @
rho2 = -0.7;

rho1 = 0.9;
rho2 = 0.0;

@ S-S Values @

kyratio=inv((gammax-beta*(1-deltak))/(beta*alp));
si=(gammax-1+deltak)*kyratio;
sc=(1-si);
etaa = (gammax-beta*(1-deltak))/gammax;
etak = (alp-1)*etaa;
etan = -etak;

@ BASIC SYSTEM MATRICES: @

scn1 = (-1  ~  0)|
       ( 0  ~  alp+gama);

scn2 = (0  ~  1  ~  0 ~ 0)|
       (alp ~  1  ~ 1 ~ 0);

smk1 = (etak  ~  1  ~  etaa ~ 0)|
       (si*gammax/(gammax-1+deltak) ~  0 ~ 0 ~ 0)|
       (0  ~  0  ~  1 ~ 0)|
       (0  ~  0  ~  0 ~ 1);

smk2 = (0  ~  1  ~  0 ~ 0)|
       (alp+si*(1-deltak)/(gammax-1+deltak)  ~  0  ~  1 ~ 0)|
       (0  ~  0    ~  rho1 ~ rho2)|
       (0  ~  0  ~  1 ~ 0);

smc1 = (0  ~  -etan)|
       (0~0)|
       (0~0)|
       (0~0);

smc2 = (0 ~ 0)|
       (-sc ~ 1-alp)|
       (0 ~ 0)|
       (0 ~ 0);

sma = 0|0|1|0;


M1 = smk1-smc1*inv(scn1)*scn2;

M2 = smk2+smc2*inv(scn1)*scn2;

W = inv(M1)*M2;

RR = inv(M1)*sma;


@ EIGENVECTOR-EIGENVALUE DECOMPOSITION OF STATE TRANSITION MATRIX:@

{vmu,vai,p,vei} = eigrg2(W);
vmu = complex(vmu,vai);
p = complex(p,vei);

@ Order eigenvalues from smallest to largest--the order of composite variables
  also change in the same way @

xx=abs(vmu)~vmu~p.';
yy=sortc(xx,1);
vmu=yy[.,2];
p=(yy[.,3:cols(yy)]).';

mu=diagrv(eye(rows(vmu)),vmu);
ps=inv(p);
/*
@ Make Sure Eigen Decomposition is correct @
ws=p*mu*ps;
if sumc(sumc(abs(ws-w))) .gt .0001;
 "Eigen decomposition of W is incorrect -- Processing stops";
  stop;
endif;
*/
"VMU = ";
VMU;


lamk = -inv(PS[4,2])*PS[4,1 3 4];

klk = W[1 3 4,1 3 4]+W[1 3 4,2]*lamk;

MKE = klk;
gm = RR[1 3 4,.];

cnk = inv(scn1)*scn2;
cnk = cnk[.,2]*lamk + cnk[.,1 3 4];

yyk = (1-alp)*cnk[2,.] + (alp ~ 1 ~ 0);
iik = 1/si*(yyk - sc*cnk[1,.]);

H = yyk|cnk[1,.]|iik|cnk[2,.];

mm = real(mke);
gm = gm*g;
H = real(H);

@ We want the variables to be:
  y, c, inv, n @

hm=zeros(4,cols(mke));
hm[1,.]=h[1,.];
hm[2,.]=h[2,.];
hm[3,.]=h[3,.];
hm[4,.]=h[4,.];

@ end of program @


*********************************************************************************************************************************
*	Wen.gss															*
*********************************************************************************************************************************

/*
   Wen.gss, YW, 01-27-2001 -- Standard KPR Model with
   -- one lag in employment decision
   -- two lags in investment decision
   -- Technology shock
   -- Government shock
*/

@ ECONOMIC PARAMETER VALUES @

alp = 0.3;       @ Capital's share @
del = 0.025;     @ Depreciation Rate for capital @
gama= 0.0;       @ > 0 singularity occurs at 0, Elasiticity of labor supply @
beta = 0.99;
sg=0.15;         @ G expenditure - output ratio @

g1 = 0*1;    @ std of tech shock @
g2 = 3;
g = (g1 ~ 0)|
    (0 ~ g2);

rhoa = 0.9;
rhog = 0.9;
RHO = (rhoa~0)|
      (0~rhog);

@ S-S Values @

kyratio=beta*alp/(1-beta*(1-del));
si=del*kyratio;
sc=(1-si-sg);
etaa = 1-beta*(1-del);
etak = (alp-1)*etaa;
etan = -etak;

@ BASIC SYSTEM MATRICES for E(c|t-2) E(n|t-2): @

scn1 = (-1  ~  0)|
       ( 0  ~  alp+gama);

scn2 = (0  ~  1)|
       (alp~  1);

scna = (0 ~ 0)|
       (1 ~ 0);

scnk = inv(scn1)*scn2;
scna = inv(scn1)*scna;

@ BASIC SYSTEM MATRICES for E[k(t+1)|t-2] E[lam(t+1)|t-2]: @

smk1 = (etak  ~  1)|
       (si/del~  0);

smk2 = (0  ~  1)|
       (alp+si*(1-del)/del  ~  0);

smc1 = (0  ~  -etan)|
       (0~0);

smc2 = (0 ~ 0)|
       (-sc ~ 1-alp);

sma1 = (-etaa ~ 0)|
       (0    ~  0);

sma2 = (0 ~ 0)|
       (1 ~ -sg);


M1 = smk1-smc1*scnk;
M2 = smk2+smc2*scnk;
R1 = sma1+smc1*scna;
R2 = sma2+smc2*scna;

W = inv(M1)*M2;
R1 = inv(M1)*R1;
R2 = inv(M1)*R2;

W = (W  ~ R1*RHO+R2)|
    (zeros(2,2)~RHO);

/* Now the state vector is E[k(t+1) lam(t+1) A(t+1) G(t+1)|t-2]
  and we need to solve for E[lam(t)|t-2] */

@ EIGENVECTOR-EIGENVALUE DECOMPOSITION OF STATE TRANSITION MATRIX:@

{vmu,vai,p,vei} = eigrg2(W);
vmu = complex(vmu,vai);
p = complex(p,vei);

@ Order eigenvalues from smallest to largest--the order of composite variables
  also change in the same way @

xx=abs(vmu)~vmu~p.';
yy=sortc(xx,1);
vmu=yy[.,2];
p=(yy[.,3:cols(yy)]).';

mu=diagrv(eye(rows(vmu)),vmu);
ps=inv(p);

@ Make Sure Eigen Decomposition is correct @
ws=p*mu*ps;
if sumc(sumc(abs(ws-w))) .gt .0001;
 "Eigen decomposition of W is incorrect -- Processing stops";
  stop;
endif;

"VMU = ";
VMU;


/* -- Step 1: Solve for E[lam(t)|t-2] in terms of E[k|t-2] E[A|t-2] E[G|t-2] */

lamk = -inv(PS[4,2])*PS[4,1 3 4];
lamk = (lamk[1] ~ rhoa^2*lamk[2] ~ rhog^2*lamk[3]);

@ Solve for E[k(t+1)|t-2] in terms of E[k|t-2] E[A|t-2] E[G|t-2] @
klk = (W[1,1]~Rhoa^2*W[1,3]~rhog^2*W[1,4])+W[1,2]*lamk;


/* -- Step 2: Solve for E[lam(t)|t-1] and E[n(t)|t-1] */
Lnn = (-sc ~ alp-1)|
      (1  ~  -alp-gama);

Lk1 = (-si/del|0);

Lk2 = (alp+si*(1-del)/del)|
      (-alp);

Lna = (1 ~ -sg)|
      (-1 ~ 0);

Ln1 = inv(Lnn)*Lk1;  @ k(t+1) @
Ln2 = inv(Lnn)*Lk2;  @ k(t)   @
Ln3 = inv(Lnn)*Lna;  @ X(t)   @
Ln3 = (rhoa*Ln3[.,1] ~ rhog*Ln3[.,2]); @ make it X(t-1) @

@ solve for E[lam(t)|t-1] E[n(t)|t-1] in terms of k(t) X(t-1) X(t-2) @
klk = klk[1] ~ 0~0~0~0 ~ klk[.,2 3];  @ make up dimension @
LNK = Ln1*klk + (Ln2~zeros(2,2)~Ln3~zeros(2,2));
nnk = LNK[2,.];

/* -- Step 3: Solve for c(t) in terms of k(t) X(t) X(t-1) X(t-2) */
cck = 1/sc*(-si/del*klk + (alp+si*(1-del)/del~1~-sg~0~0~0~0)+(1-alp)*nnk );


/* -- Step 4: Solve for y(t) */
yyk = (alp~1~0~0~0~0~0) + (1-alp)*nnk;

iik = 1/del*klk - (1-del)/del*(1 ~ zeros(1,6));
iik1 = 1/si*(yyk - sc*cck - sg*(0~0~1~zeros(1,4)));

/* -- Step 5: arrange the state vector k(t) X(t) X(t-1) X(t-2) */
mke = (klk)|
     (zeros(2,1)~RHO~zeros(2,4))|
     (zeros(2,1)~eye(2)~zeros(2,4))|
     (zeros(2,3)~eye(2)~zeros(2,2));

gm = (zeros(1,2)|eye(2)|zeros(4,2));

H = yyk|cck|iik|nnk;

mm = real(mke);
gm = gm*g;
H = real(H);

@ We want the variables to be:
  y, c, inv, n @

hm=zeros(4,cols(mke));
hm[1,.]=h[1,.];
hm[2,.]=h[2,.];
hm[3,.]=h[3,.];
hm[4,.]=h[4,.];

@ end of program @


*********************************************************************************************************************************
*	Granger.gss														*
*********************************************************************************************************************************

/*
   Granger.gss, YW, 02-10-2001 -- Benhabib-Wen Model with
   -- one lag in employment decision
   -- two lags in investment decision
   -- lagged Capacity utilization cost
   -- Technology shock
   -- Government shock
*/

@ ECONOMIC PARAMETER VALUES @

eta1 = 0.15;     @ Externality @
alp = 0.3;       @ Capital's share @
del = 0.025;     @ Depreciation Rate for capital @
gama= 0.0;       @ > 0 singularity occurs at 0, Elasiticity of labor supply @
beta = 0.99;
sg=0.15;         @ G expenditure - output ratio @

g1 = 0.2;    @ std of tech shock @
g2 = 1;
g = (0*g1 ~ 0)|
    (0 ~ g2);

rhoa = 0.9;
rhog = 0.9;
RHO = (rhoa ~ 0)|
      (0 ~ rhog);

@ S-S Values @

kyratio=beta*alp/(1-beta*(1-del));
the = (1-beta*(1-del))/(beta^2*del);
si=del*kyratio;
sc=(1-si-sg);
etaa = 1-beta*(1-del);
etak = (alp*(1+eta1)-1)*etaa;
etan = (1-alp)*(1+eta1)*etaa;
etae = alp*(1+eta1)*etaa;

@ BASIC SYSTEM MATRICES for E(c|t-2) E(n|t-2) @

scn1 = (-1  ~  0)|
       ( 0  ~  1+gama-(1-alp)*(1+eta1));

scn2 = (0 ~ 1 ~ 0 ~ 0 ~ 0 ~ 0)|
       (alp*(1+eta1)~  1 ~ alp*(1+eta1) ~ 0 ~ 1 ~ 0);

scnk = inv(scn1)*scn2;

@ BASIC SYSTEM MATRICES for k(t+1) lam(t+1) e(t+1) e(t) X(t+1) @

smk1 = (1 ~ 1 ~ 0 ~ 0 ~ 0 ~ 0)|
       (etak ~ 1 ~ etae ~ 0 ~ etaa ~ 0)|
       (si/del~  0 ~ 0 ~ 0 ~ 0 ~ 0)|
       (0 ~ 0 ~ 0 ~ 1 ~ 0 ~ 0)|
       (0 ~ 0 ~ 0 ~ 0 ~ 1 ~ 0)|
       (0 ~ 0 ~ 0 ~ 0 ~ 0 ~ 1);

smk2 = (alp*(1+eta1) ~ 1 ~ alp*(1+eta1)-the ~ 0 ~ 1 ~ 0)|
       (0 ~ 1 ~ beta*del*the ~ 0 ~ 0 ~ 0)|
       (alp*(1+eta1)+si*(1-del)/del ~ 0 ~ alp*(1+eta1) ~ -si*the ~ 1 ~ -sg)|
       (0 ~ 0 ~ 1 ~ 0 ~ 0 ~ 0)|
       (0 ~ 0 ~ 0 ~ 0 ~ rhoa ~ 0)|
       (0 ~ 0 ~ 0 ~ 0 ~ 0 ~ rhog);

smc1 = (0  ~  0)|
       (0  ~ -etan)|
       (zeros(4,2));

smc2 = (0 ~ (1-alp)*(1+eta1))|
       (0 ~ 0)|
       (-sc ~ (1-alp)*(1+eta1))|
       (zeros(3,2));

M1 = smk1-smc1*scnk;
M2 = smk2+smc2*scnk;

W = inv(M1)*M2;

/* Now the state vector is E[k(t+1) lam(t+1) A(t+1) G(t+1)|t-2]
  and we need to solve for E[lam(t)|t-2] */

@ EIGENVECTOR-EIGENVALUE DECOMPOSITION OF STATE TRANSITION MATRIX:@

{vmu,vai,p,vei} = eigrg2(W);
vmu = complex(vmu,vai);
p = complex(p,vei);

@ Order eigenvalues from smallest to largest--the order of composite variables
  also change in the same way @

xx=abs(vmu)~vmu~p.';
yy=sortc(xx,1);
vmu=yy[.,2];
p=(yy[.,3:cols(yy)]).';

mu=diagrv(eye(rows(vmu)),vmu);
ps=inv(p);

@ Make Sure Eigen Decomposition is correct @
ws=p*mu*ps;
if sumc(sumc(abs(ws-w))) .gt .0001;
 "Eigen decomposition of W is incorrect -- Processing stops";
  stop;
endif;

"VMU = ";
VMU;


/* -- Step 1: Solve for E[lam(t),e(t)] in term E[k|t-2] E[e(t-1|t-2] E[X|t-2] */

lamk = -inv(PS[5 6,2 3])*PS[5 6,1 4 5 6];
lamk = (lamk[.,1 2] ~ rhoa^2*lamk[.,3] ~ rhog^2*lamk[.,4]);

@ Solve for E[k(t+1)|t-2] in terms of E[k|t-2] E[e(t-1)|t-2] E[X|t-2] @
klk = (W[1,1 4]~rhoa^2*W[1,5]~rhog^2*W[1,6])+W[1,2 3]*lamk;


/* -- Step 2: Solve for E[lam(t), n(t), e(t)|t-1] */

klk = klk[.,1 2] ~ 0~0~0~0 ~ klk[.,3 4];  @ make up dimension @

lne = (-1 ~ 1+gama-(1-alp)*(1+eta1) ~ -alp*(1+eta1))|
      (1 ~ (1-alp)*(1+eta1) ~ alp*(1+eta1)-the-lamk[1,2])|
      (-sc ~ -(1-alp)*(1+eta1) ~ -alp*(1+eta1));

lnek1 = (0)|(1+lamk[1,1])|(-si/del);
lnek1 = lnek1*klk;

lnek2 = (alp*(1+eta1) ~ 0 ~ rhoa ~ 0)|
        (-alp*(1+eta1) ~ 0 ~ -rhoa+lamk[1,3] ~ lamk[1,4])|
        (alp*(1+eta1)+si*(1-del)/del ~ -si*the ~ rhoa ~ -sg*rhog);
lnek2 = (lnek2[.,1 2] ~ zeros(3,2) ~ lnek2[.,3 4] ~ zeros(3,2));

LNEK = inv(lne)*(lnek1 + lnek2);
nnk = LNEK[2,.];
eek = LNEk[3,.];

/* -- Step 3: Solve for c(t) in terms of k(t) e(t-1) X(t) X(t-1) X(t-2) */
cck=-si/del*klk+(alp*(1+eta1)+si*(1-del)/del~-si*the~1~-sg~0~0~0~0);
cck = (1/sc)*(cck + (0 ~ (1-alp)*(1+eta1) ~ alp*(1+eta1))*Lnek);


/* -- Step 4: Solve for y(t) */
yyk= (alp*(1+eta1)~0~1~0~0~0~0~0) + ((1-alp)*(1+eta1)~alp*(1+eta1))*Lnek[2 3,.];

iik = 1/del*klk - (1-del)/del*(1 ~ -del*the/(1-del) ~ zeros(1,6));
iik1 = 1/si*(yyk - sc*cck - sg*(0~0~0~1~zeros(1,4)));

/* -- Step 5: arrange the state vector k(t) e(t-1) X(t) X(t-1) X(t-2) */
mke = (klk)|
      (eek)|
      (zeros(2,2)~RHO~zeros(2,4))|
      (zeros(2,2)~eye(2)~zeros(2,4))|
      (zeros(2,4)~eye(2)~zeros(2,2));

gm = (zeros(2,2)|eye(2)|zeros(4,2));

H = yyk|cck|iik|nnk;

mm = real(mke);
gm = gm*g;
H = real(H);

@ We want the variables to be:
  y, c, inv, n @

hm=zeros(4,cols(mke));
hm[1,.]=h[1,.];
hm[2,.]=h[2,.];
hm[3,.]=h[3,.];
hm[4,.]=h[4,.];

@ end of program @


*********************************************************************************************************************************
*********************************************************************************************************************************
*********************************************************************************************************************************