@GAUSS PROGRAM@

library pgraph;

@NOTE: VALUES OF qc, qy, qx DO NOT MATTER FOR ROOTS@

al1=0.1; a1=0.00; al2=0.25; a2=0.00; a0=0.00; al0=1-al1-a1-al2-a2-a0;qc=1;           @values in Working Paper@
be1=0.25; b1=0.09; be2=0.20; b2=0.00; b0=0.00; be0=1-be1-b1-be2-b2-b0;qx=1;
ga1=0.23; c1=0.00; ga2=0.10; c2=0.1; c0=0.00; ga0=1-ga1-c1-ga2-c2-c0;qy=1;
r=.036;g=.000025;s=1;nu=100;
ewc=.018;ewx=.018;ewy=.018;swx=.018;swy=.018;
lamc=0.95;lamx=.95;lamy=.95;  

             
                 
al1=0.15;   a1=0.07; al2=0.20;  a2=0.00; a0=0.00; al0=1-al1-a1-al2-a2-a0;     qc=1;           @new values @
be1=0.22; b1=0.07; be2=0.21; b2=0.00; b0=0.00; be0=1-be1-b1-be2-b2-b0;  qx=1;
ga1=0.26; c1=0.00; ga2=0.15; c2=0.08; c0=0.00; ga0=1-ga1-c1-ga2-c2-c0;    qy=1;
r=.036;g=.025;s=1;nu=5;
ewc=.0039;ewx=.0039;ewy=.0039;swx=.0039;swy=.0039;
lamc=0.95;lamx=.95;lamy=.95;                   

/*


al1=0.32; a1=0.100; al2=0.20; a2=0.00; a0=0.00; al0=1-al1-a1-al2-a2-a0;qc=1;           @values in Perli Paper@
be1=0.24; b1=0.0; be2=0.20; b2=0.20; b0=0.00; be0=1-be1-b1-be2-b2-b0;qx=1;
ga1=0.22; c1=0.0; ga2=0.05; c2=0.0; c0=0.00; ga0=1-ga1-c1-ga2-c2-c0;qy=1;
r=.036;g=.025;s=1;nu=100;
ewc=.018;ewx=.018;ewy=.018;swx=.018;swy=.018;
lamc=0.95;lamx=.95;lamy=.95;                                    


								@no externalities and therefore no indeterminacy@

al1=0.15; a1=0.00; al2=0.20; a2=0.00; a0=0.00; al0=1-al1-a1-al2-a2-a0;qc=1;	
be1=0.22; b1=0.0; be2=0.21; b2=0.00; b0=0.00; be0=1-be1-b1-be2-b2-b0;qx=1;
ga1=0.26; c1=0.00; ga2=0.15; c2=0; c0=0.00; ga0=1-ga1-c1-ga2-c2-c0;qy=1;
r=.036;g=.025;s=1;nu=10;
ewc=.018;ewx=.018;ewy=.018;swx=.018;swy=.018;
lamc=0.95;lamx=.95;lamy=.95;       



al1=0.11; a1=0.00; al2=0.25; a2=0.00; a0=0.00; al0=1-al1-a1-al2-a2-a0;qc=1;		@other useful values@
be1=0.28; b1=0.07; be2=0.20; b2=0.00; b0=0.00; be0=1-be1-b1-be2-b2-b0;qx=1;
ga1=0.23; c1=0.00; ga2=0.10; c2=0.07; c0=0.00; ga0=1-ga1-c1-ga2-c2-c0;qy=1;
r=.036;g=.025;s=1;nu=1;
ewc=.018;ewx=.018;ewy=.018;swx=.018;swy=.018;
lamc=0.95;lamx=.95;lamy=.95;

*/


@ nu is the labor supply elasticity. nu=0 implies fixed labor supply, while
nu=infinity is utility linear in labor or leisure. It is equal to 1/v, where
the utility function is U(c,L)=(1/(1-s))*c^(1-s)-(1/(1+v))*L^(1+v). The
Cobb-Douglas exponents alpha and a, (al0, al1, al2, a0, a1, a2) are for the
consumptiongood, where the a's are the external effects, The beta's
and b's (be0, be1, be2, b0, b1, b0) are for the investment good x,and the
gamma's and c's, (ga0, ga1, ga2, c0, c1, c2) are for the investment good y.@

    @vc matrix for shocks:  @


omega=eye(6); 														@independent shocks@
omega=eye(2)~zeros(2,4)|zeros(4,2)~ones(4,4);				@perfectly  correlated technology / sunspot shocks@	
omega=eye(4)~zeros(4,2)|zeros(2,4)~ones(2,2);				@perfectly  correlated technology shocks@
omega=eye(3)~zeros(3,3)|zeros(3,3)~ones(3,3);

@STEADY STATE PRICES@

sp0=be0+b0~be2+b2|ga0+c0~-ga1-c1-ga0-c0;

lcs0=ln(r)-ln(qx*be1*((be0/be1)^(be0+b0))*((be2/be1)^(be2+b2)));

lcs1=ln(r)-ln(qy*ga2*((ga1/ga2)^(ga1+c1))*((ga0/ga2)^(ga0+c0)));

lcs2= ln((qc/qx)*(al1/be1)*((al2/al1)^(al2+a2))*((al0/al1)^(al0+a0))
*((be2/be1)^(-be2-b2))*((be0/be1)^(-be0-b0)));

lcs3= ln((qc/qy)*(al1/ga1)*((al2/al1)^(al2+a2))*((al0/al1)^(al0+a0))
*((ga2/ga1)^(-ga2-c2))*((ga0/ga1)^(-ga0-c0)));

sp1=al0+a0-be0-b0~al2+a2-be2-b2|al0+a0-ga0-c0~al2+a2-ga2-c2;

lsom=inv(sp0)*(lcs0|lcs1);

lsp=(lcs2|lcs3)+(sp1)*lsom;

om=exp(lsom);

ps=exp(lsp);

wx=r*ps[1];

wy=r*ps[2];

IM=1~0|0~1;

lom=inv(sp1)*(lsp-(lcs2|lcs3));

lw=(-lcs0+ln(r)|-lcs1+ln(r))+(sp0*inv(sp1)+IM)*lsp-sp0*inv(sp1)*(lcs2|lcs3);

@ INPUT COEFFICIENTS @

a01=(1/qx)*(((be0/be1)*om[1])^(be1+b1))*(((be0/be2)*(om[1]/om[2]))^(be2+b2));
a11=(1/qx)*(((1/om[1])*(be1/be0))^(be0+b0))*(((be1/be2)*(1/om[2]))^(be2+b2));
a21=(1/qx)*(((be2/be1)*om[2])^(be1+b1))*(((be2/be0)*(om[2]/om[1]))^(be0+b0));

a02=(1/qy)*(((ga0/ga1)*om[1])^(ga1+c1))*(((ga0/ga2)*(om[1]/om[2]))^(ga2+c2));
a12=(1/qy)*(((1/om[2])*(ga1/ga2))^(ga2+c2))*(((ga1/ga0)*(1/om[1]))^(ga0+c0));
a22=(1/qy)*(((ga2/ga0)*(om[2]/om[1]))^(ga0+c0))*(((ga2/ga1)*om[2])^(ga1+c1));

a00=(1/qc)*(((al0/al1)*om[1])^(al1+a1))*(((al0/al2)*(om[1]/om[2]))^(al2+a2));
a10=(1/qc)*(((1/om[1])*(al1/al0))^(al0+a0))*(((al1/al2)*(1/om[2]))^(al2+a2));
a20=(1/qc)*(((al2/al0)*(om[2]/om[1]))^(al0+a0))*(((al2/al1)*om[2])^(al1+a1));

@ STEADY STATE QUANTITIES @

A=a11~a12|a21~a22;

cs=((a00+g*(a01~a02)*(inv(IM-g*A))*(a10|a20))*((a00/al0)^nu))^(-1/(nu*s+1));

kxys=cs*(inv(IM-g*A)*(a10|a20));
ys=g*kxys[2];
xs=g*kxys[1];

@ ELASTICITIES: @

lc=a00*cs; kxc=a10*cs; kyc=a20*cs;
lx=a01*xs; kxx=a11*xs; kyx=a21*xs;
ly=a02*ys; kxy=a12*ys; kyy=a22*ys;

kx=kxc+kxx+kxy;
ky=kyc+kyx+kyy;
l=lc+lx+ly;

x1=(al1+a1)*(1-s)*nu;
x2=(al2+a2)*(1-s)*nu;
x0=((al0+a0)*(1-s)*nu)-nu;

F=((lc/l)+(s*nu))~lx/l~ly/l|kxc/kx~kxx/kx~kxy/kx|kyc/ky~kyx/ky~kyy/ky;
FI=inv(F);
SI=inv(sp1);

khom12h=FI*(1-x0-x1|1|0);
khom10h=FI*(-1+x0|0|0);

khkxh=FI*(0|1|0);
khkyh=FI*(0|0|1);

khpxh=FI*((-SI[1,1]+SI[2,1])*(1-x0)-x1*SI[2,1]|SI[2,1]|0);
khpyh=FI*((-SI[1,2]+SI[2,2])*(1-x0)-x1*SI[2,2]|SI[2,2]|0);

chpx=khpxh[1]+(al0+a0)*(SI[1,1]-SI[2,1])-(al1+a1)*SI[2,1];
chpy=khpyh[1]+(al0+a0)*(SI[1,2]-SI[2,2])-(al1+a1)*SI[2,2];
chkx=khkxh[1];
chky=khkyh[1];

xhpx=khpxh[2]+(be0+b0)*(SI[1,1]-SI[2,1])-(be1+b1)*SI[2,1];
xhpy=khpyh[2]+(be0+b0)*(SI[1,2]-SI[2,2])-(be1+b1)*SI[2,2];
xhkx=khkxh[2];
xhky=khkyh[2];

yhpx=khpxh[3]+(ga0+c0)*(SI[1,1]-SI[2,1])-(ga1+c1)*SI[2,1];
yhpy=khpyh[3]+(ga0+c0)*(SI[1,2]-SI[2,2])-(ga1+c1)*SI[2,2];
yhkx=khkxh[3];
yhky=khkyh[3];

@ JACOBIAN @
px=ps[1];
py=ps[2];

J11=(xhkx*xs/kx~xhky*xs/ky|yhkx*ys/kx~yhky*ys/ky)-(g*IM);
J12=(xhpx*xs/px~xhpy*xs/py|yhpx*ys/px~yhpy*ys/py);

D=IM-s*(px/cs~0|0~py/cs)*(chpx*cs/px~chpy*cs/py|chpx*cs/px~chpy*cs/py);
DI=inv(D);

J21=DI*(s*(px/cs~0|0~py/cs))*(chkx*cs/kx~chky*cs/ky|chkx*cs/kx~chky*cs/ky)*
((xhkx*xs/kx~xhky*xs/ky|yhkx*ys/kx~yhky*ys/ky)-(g*IM));

WW=sp0*inv(sp1)+IM;

WR=WW[1,1]*r~WW[1,2]*wx/py|WW[2,1]*wy/px~WW[2,2]*r;

J22=DI*(s*(px/cs~0|0~py/cs))*(chkx*cs/kx~chky*cs/ky|chkx*cs/kx~chky*cs/ky)*
(xhpx*xs/px~xhpy*xs/py|yhpx*ys/px~yhpy*ys/py)+DI*(-WR+r*IM);
MMM=(s*(px/cs~0|0~py/cs))*(chkx*cs/kx~chky*cs/ky|chkx*cs/kx~chky*cs/ky);

J=J11~J12|J21~J22;

@ROOTS FOR CONTINOUS TIME@

roots=eig(J);
print "roots for continous time case";roots';


@ SOLVING DISCRETE TIME CASE WITH SUNSPOT SHOCKS@

IM3=1~0~0|0~1~0|0~0~1;
z2=zeros(2,2);
z23=zeros(2,3);
z3=zeros(3,3);

zxbar=1;zybar=1;zcbar=1;



wwop=WW-IM;
RD11=((xhkx*xs/kx)+(1-g))*kx/(xs+(1-g)*kx)~xhky*xs/(xs+(1-g)*kx)
|yhkx*ys/(ys+(1-g)*ky)~(ky/(ys+(1-g)*ky))*((yhky*ys/ky)+(1-g));
RD12=xhpx*xs/(xs+(1-g)*kx)~xhpy*xs/(xs+(1-g)*kx)|yhpx*ys/(ys+(1-g)*ky)~
yhpy*ys/(ys+(1-g)*ky);
RD21=-s*chkx~-s*chky|-s*chkx~-s*chky;
RD22=-s*chpx+1~-s*chpy|-s*chpx~-s*chpy+1;
RD13=(xs/kx)*(-xhpx-xhpy)~(xs/kx)*(1+xhpx)~(xs/kx)*(xhpy)|
(ys/ky)*(-yhpx-yhpy)~(ys/ky)*(yhpx)~(ys/ky)*(1+yhpy);
RD23=-s*(1-chpx-chpy)~-s*(chpx)~-s*(chpy)|
     -s*(1-chpx-chpy)~-s*(chpx)~-s*(chpy);
RD31=0~0|0~0|0~0;
RD32=RD31;
RD33=lamc~0~0|0~lamx~0|0~0~lamy;

RD=RD11~RD12~RD13|RD21~RD22~RD23|RD31~RD32~RD33;

QD11=1~0|0~1;
QD12=0~0|0~0;
QD13=0~0~0|0~0~0;

QD21=-s*chkx~-s*chky|-s*chkx~-s*chky;

QD22=(1/(1+r-g))*r*wwop[1,1]-s*chpx+1~(1/(1+r-g))*r*wwop[1,2]-s*chpy|
(1/(1+r-g))*r*wwop[2,1]-s*chpx~(1/(1+r-g))*r*wwop[2,2]-s*chpy+1;

QD23=(1/(1+r-g))*r*(-wwop[1,1]-wwop[1,2])-s*(1-chpx-chpy)
~(1/(1+r-g))*(1+wwop[1,1])-s*chpx
~(1/(1+r-g))*r*wwop[1,2]-s*chpy|
 (1/(1+r-g))*r*(-wwop[2,1]-wwop[2,2])-s*(1-chpx-chpy)
~(1/(1+r-g))*r*(wwop[2,1])-s*chpx
~(1/(1+r-g))*r*(wwop[2,2]+1)-s*chpy;

QD31=RD31;
QD32=RD31;

QD33=1~0~0|0~1~0|0~0~1;

QD=QD11~QD12~QD13|QD21~QD22~QD23|QD31~QD32~QD33;

JD=inv(QD)*RD;

Print"";
Format 4,3;
print "Roots for discrete time: ";eig(JD)';
{va,ve}=eigv(JD);
modva=sqrt(real(va)^2+imag(va)^2);
vv=inv(ve);
interv={1, 1e12};
inx=indexcat(modva,interv);

if rows(inx)==1;

jjd=jd[1:3,.]|jd[5:7,.];
vdy=-(1/vv[inx[1],4])*(vv[inx[1],1:3]~vv[inx[1],5:7]);

njjd1=jjd[.,1:3]~jjd[.,5:7];
njjd2=((jd[1:3,4]|jd[5:7,4])*vdy);
nj=njjd1+njjd2;






elseif rows(inx)== 2;

jjd=jd[1:2,.]|jd[5:7,.];
vvsub=submat(vv,inx,0);
vvr=-inv(vvsub[.,3:4])*(vvsub[.,1:2]~vvsub[.,5:7]);


njjd1=jjd[.,1:2]~jjd[.,5:7];
njjd2=((jd[1:2,3:4]|jd[5:7,3:4])*vvr);
nj=njjd1+njjd2;

else;
print "more than two roots with modulus > 1";
end;

endif;


yt=0|0|0|0|0|0;
y=yt|0|0|0|0|0|0;


IM3=1~0~0|0~1~0|0~0~1;
z2=zeros(2,2);
z23=zeros(2,3);
z3=zeros(3,3);

BM=0~0~0~-1~0~0~1~0~0| 1~0~0~-1~0~0~0~0~0|0~0~0~0~-1~0~0~1~0|0~1~0~0~-1~0~0~0~0|
0~0~0~0~0~-1~0~0~1|0~0~1~0~0~-1~0~0~0|
0~0~0~kxc/kx~kxx/kx~kxy/kx~0~0~0|
kyc/ky~kyx/ky~kyy/ky~0~0~0~0~0~0|
-x2~0~0~-x1~0~0~lc/l-x0~lx/l~ly/l;

NN=inv(sp1)~z2~z2~z23|z2~inv(sp1)~z2~z23|z2~z2~inv(sp1)~z23|z3~z3~IM3;

SS=inv(BM)*NN;

gnp=cs+px*xs+py*ys;
cgnp=cs/gnp;
xgnp=px*xs/gnp;
ygnp=py*ys/gnp;


@ANALYTIC MOMENTS@

    @vc matrix for shocks -I put it with the parameters at the  beginning@

if rows(inx)==1;



@omega=1~0~0~0~0~0|0~1~0~0~0~0|0~0~1~0~0~0|0~0~0~1~0~0|0~0~0~0~1~0|0~0~0~0~0~1;@

                    @coefficients for shocks@
q2=eye(6);
q2[1,1]=0;
q2[2,2]=0;
q2[3,3]=swx;
q2[4,4]=ewc;
q2[5,5]=ewx;
q2[6,6]=ewy;

                     @vc matrix for shocks with coefficients@
varqom0=q2*omega*q2';

QDI=inv(QD);
qqd=QDI[1:3,.]|QDI[5:7,.];
qqdv=qqd[.,1:3]~qqd[.,5:7];

VQD=vv*QDI;
ssyv=-(VQD[inx[1],4]^(-1))*(VQD[inx[1],1:3]~VQD[inx[1],5:7]);
qqdd=qqdv+((QDI[1:3,4]|QDI[5:7,4])*ssyv);


varqom=qqdd*varqom0*qqdd';

                @vc for kx,ky,px,zc,zx,zy@
kk=nj.*.nj;
cvos=inv(eye(36)-kk)*vec(varqom);
cv0=reshape(cvos,6,6);
 @add py to vc matrix. The sequence is critical to understand what follows@
                   @vc for py,kx,ky,px,zc,zx,zy@

cv0py=(vdy|eye(6))*cv0*(vdy'~eye(6));

       @add inputs,kyc,kyx,kyy,kxc,kxx,kxy,Lc,Lx,Ly to vc matrix@

jj0=
0~0~0~1~0~0~0|1~0~0~0~0~0~0|
0~0~0~1~0~0~0|1~0~0~0~0~0~0|
0~0~0~1~0~0~0|1~0~0~0~0~0~0|
0~1~0~0~0~0~0|0~0~1~0~0~0~0|0~0~0~0~0~0~0;
jj1=SS*jj0;
jji=eye(7)|jj1;

cvxy=jji*cv0py*jji';

@cvxy is vc matrix is without outputs.  We will add them
Note sequence: py,kx,ky,px,zc,zx,zy,kyc,kyx,kyy,kxc,kxx,kxy,Lc,Lx,Ly.
Below is vc matrix including output growth of c, x, y with shoks and prices@
outp=
0~0~0~0~1~0~0~al2+a2~0~0~al1+a1~0~0~al0+a0~0~0|
0~0~0~1~0~1~0~0~be2+b2~0~0~be1+b1~0~0~be0+b0~0|
1~0~0~0~0~0~1~0~0~ga2+c2~0~0~ga1+c1~0~0~ga0+c0;

ncvxy=(outp|eye(16))*cvxy*(outp|eye(16))';


                @now add gnp @
@Note new sequence: c,px*x,py*y,py,kx,ky,px,zc,zx,zy,kyc,kyx,kyy,kxc,kxx,kxy,
Lc,Lx,Ly.@

vgnp=(cgnp~xgnp~ygnp~0~0~0~0~0~0~0~0~0~0~0~0~0~0~0~0);


ougn=vgnp|eye(19);
cvgnp=ougn*ncvxy*ougn';

                                       @vc for gnp,c,px*x,
                                    py*y,py,kx,ky,px.zc,zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly@
    @Now add total investment, to the end of the string this time@

invsum=px*xs+py*ys;
sixs=px*xs/invsum;
siys=py*ys/invsum;
vsi=(0~0~sixs~siys);
zb=zeros(1,16);
vsiv=eye(20)|(vsi~zb);
cvgnpv=vsiv*cvgnp*vsiv';

                                                @vc for gnp,c,px*x,
                                    py*y,py,kx,ky,px,zc.zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly,(px*x+py*y)@

                  @Finally add total labor@

tl0=zb~0~lc/l~lx/l~ly/l~0;
tl=eye(21)|tl0;
cvgnpz=tl*cvgnpv*tl' ;

                                                     @cv for gnp,c,px*x,
                                    py*y,py,kx,ky,px.zc.zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly,(px*x+py*y),l@
cvgnpz=real(cvgnpz);

cv=corrvc(cvgnpz);                            @correlation matrix@

@ AUTOCORROLATION COMPUTATION AT ONE LAG. FOR LAG OF ORDER t, SIMPLY MULTIPLY THE
MATRIX nj IN THE EQUATION JUST BELOW t TIMES BY ITSELF (note raising a matrix to a power does 
not work because gauss does elememt by element multiplication when a matrix is raised to a power . @

nj=real(nj);
cv0=real(cv0);
arcv0=(nj)*cv0;
arcv0py=(vdy|eye(6))*arcv0*(vdy'~eye(6));
arcvxy=jji*arcv0py*jji';
arncvxy=(outp|eye(16))*arcvxy*(outp|eye(16))';
arcvgnp=ougn*arncvxy*ougn';
arcvgnpv=vsiv*arcvgnp*vsiv';
arcvgnpz=tl*arcvgnpv*tl' ;

arcv=sqrt((inv(diagrv(eye(22),diag(cvgnpz)))))*arcvgnpz
*sqrt((inv(diagrv(eye(22),diag(cvgnpz)))));


                        @standard deviations@

stdcs=sqrt(cvgnpz[2,2]);
stdx=sqrt(cvgnpz[3,3]);
stdy=sqrt(cvgnpz[4,4]);
stdlc=sqrt(cvgnpz[18,18]);
stdkx=sqrt(cvgnpz[6,6]);
stdky=sqrt(cvgnpz[7,7]);
stdzc=sqrt(cvgnpz[9,9]);
stdzx=sqrt(cvgnpz[10,10]);
stdgnp=sqrt(cvgnpz[1,1]);
sdti=sqrt(cvgnpz[21,21]);
stdl=sqrt(cvgnpz[22,22]);
stdpy=sqrt(cvgnpz[5,5]);
stdpx=sqrt(cvgnpz[8,8]);

                @relative std deviations@

rstdx=stdx/stdgnp;
rstdc=stdcs/stdgnp;
rstdlc=stdlc/stdgnp;
rstdkx=stdkx/stdgnp;
rstdzc=stdzc/stdgnp;
rstdy=stdy/stdgnp;
rstdi=sdti/stdgnp;
rstdl=stdl/stdgnp;
cv2=cv[1,2];cv3=cv[1,3];cv21=cv[1,21];cv22=cv[1,22];cv4=cv[1,4];

            @Printing results@

print"
--------------------------- Standard Deviations ------------------------------";
print "        Output          C               X                  I                L               Y";
format/rd 12,4;
print "";
print 1~rstdc~rstdx~rstdi~rstdl~rstdy;
print "           (-)            (.73)            ( -)              (3.20)         (1.16)           (-)
";
print"
------------------------ Correlations with Output ----------------------------";
print "      Output            C               X                 I                  L               Y"
;
format 12,4;
print "";
print 1~cv[1,2]~cv[1,3]~cv[1,21]~cv[1,22]~cv[1,4];
print "         (-)              (.82)             (-)             (.90)            ( .86 )             (-)
";
print"
--------------------------- AR(1) Coefficients -------------------------------";
print "Output       C             X              I              L                Y";

format 4,2;
print  arcv[1,1]~arcv[2,2]~arcv[3,3]~arcv[21,21]~arcv[22,22]~arcv[4,4];

print "(.90)         (.84)          (-)            (.76)       ( .90 )            (-)";

                                                     @cv for gnp,c,px*x,
                                    py*y,py,kx,ky,px.zc.zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly,(px*x+py*y),l@



elseif rows(inx)==2;




omega=1~0~0~0~0|0~1~0~0~0|0~0~1~0~0|0~0~0~1~0|0~0~0~0~1;

                    @coefficients for shocks@
q2=eye(5);
q2[1,1]=0;
q2[2,2]=0;
q2[3,3]=ewc;
q2[4,4]=ewx;
q2[5,5]=ewy;

                     @vc matrix for shocks with coefficients@
varqom0=q2*omega*q2';
QDI=inv(QD);
VQD=vv*QDI;
vqdsub=submat(vqd,inx,0);

ssyx=-inv(vqdsub[.,3:4])*(vqdsub[.,1:2]~vqdsub[.,5:7]);
qqd=QDI[1:2,.]|QDI[5:7,.];
qqdv=qqd[.,1:2]~qqd[.,5:7];
qqdd=qqdv+(QDI[1:2,3:4]|QDI[5:7,3:4])*ssyx;
varqom=qqdd*varqom0*qqdd';


                @vc for kx,ky,zc,zx,zy@
kk=nj.*.nj;

sqs=inv(eye(25)-kk);

cvss0=sqs*vec(varqom);


cv0=reshape(cvss0,5,5);

 @add px,py to vc matrix. The sequence is critical to understand what follows@
                   @vc for px,py,kx,ky,zc,zx,zy@

cv0py=(vvr|eye(5))*cv0*(vvr'~eye(5));

        @NOTE; SEQUENCE IS SLIGHTLY DIFFERENT BECAUSE OF PX@
       @add inputs,kyc,kyx,kyy,kxc,kxx,kxy,Lc,Lx,Ly to vc matrix@

jj0=
1~0~0~0~0~0~0|0~1~0~0~0~0~0|
1~0~0~0~0~0~0|0~1~0~0~0~0~0|
1~0~0~0~0~0~0|0~1~0~0~0~0~0|
0~0~1~0~0~0~0|0~0~0~1~0~0~0|0~0~0~0~0~0~0;
jj1=SS*jj0;
jji=eye(7)|jj1;

cvxy=jji*cv0py*jji';

@cvxy is vc matrix is without outputs.  We will add them
Note sequence: px,py,kx,ky,zc,zx,zy,kyc,kyx,kyy,kxc,kxx,kxy,Lc,Lx,Ly.
Below is vc matrix including output growth of c, x, y with shoks and prices@
outp=
0~0~0~0~1~0~0~al2+a2~0~0~al1+a1~0~0~al0+a0~0~0|
1~0~0~0~0~1~0~0~be2+b2~0~0~be1+b1~0~0~be0+b0~0|
0~1~0~0~0~0~1~0~0~ga2+c2~0~0~ga1+c1~0~0~ga0+c0;

ncvxy=(outp|eye(16))*cvxy*(outp|eye(16))';

                @now add gnp @
@Note new sequence: c,px*x,py*y,px,py,kx,ky,zc,zx,zy,kyc,kyx,kyy,kxc,kxx,kxy,
Lc,Lx,Ly.@

vgnp=(cgnp~xgnp~ygnp~0~0~0~0~0~0~0~0~0~0~0~0~0~0~0~0);

ougn=vgnp|eye(19);
cvgnp=ougn*ncvxy*ougn';                       @vc for gnp,c,px*x,
                                    py*y,px,py,kx,ky,zc,zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly@
    @Now add total investment, to the end of the string this time@

invsum=px*xs+py*ys;
sixs=px*xs/invsum;
siys=py*ys/invsum;
vsi=(0~0~sixs~siys);
zb=zeros(1,16);
vsiv=eye(20)|(vsi~zb);
cvgnpv=vsiv*cvgnp*vsiv';           @vc for gnp,c,px*x,
                                    py*y,px,py,kx,ky,zc.zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly,(px*x+py*y)@

                  @Finally add total labor@

tl0=zb~0~lc/l~lx/l~ly/l~0;
tl=eye(21)|tl0;
cvgnpz=tl*cvgnpv*tl' ;                   @cv for gnp,c,px*x,
                                    py*y,px,py,kx,ky,zc.zx,zy,kyc,kyx,kyy,
                                             kxc,kxx,kxc,lc,lx,ly,(px*x+py*y),l@
cvgnpz=real(cvgnpz);

cv=corrvc(cvgnpz);                            @correlation matrix@

                @AUTOCORROLATION COMPUTATION@

@ AUTOCORROLATION COMPUTATION AT ONE LAG. FOR LAG OF ORDER t, SIMPLY MULTIPLY THE
MATRIX nj IN THE EQUATION JUST BELOW t TIMES BY ITSELF (note raising a matrix to a power does 
not work because gauss does elememt by element multiplication when a matrix is raised to a power . @

nj=real(nj);
cv0=real(cv0);
arcv0=nj*cv0;
arcv0py=(vvr|eye(5))*arcv0*(vvr'~eye(5));
arcvxy=jji*arcv0py*jji';
arncvxy=(outp|eye(16))*arcvxy*(outp|eye(16))';
arcvgnp=ougn*arncvxy*ougn';
arcvgnpv=vsiv*arcvgnp*vsiv';
arcvgnpz=tl*arcvgnpv*tl' ;

arcv=sqrt((inv(diagrv(eye(22),diag(cvgnpz)))))*arcvgnpz
*sqrt((inv(diagrv(eye(22),diag(cvgnpz)))));


             @CAUTION:ORDERING IS DIFFERENT THAN INDET. CASE@
                        @standard deviations@

stdcs=sqrt(cvgnpz[2,2]);
stdx=sqrt(cvgnpz[3,3]);
stdy=sqrt(cvgnpz[4,4]);
stdlc=sqrt(cvgnpz[18,18]);
stdkx=sqrt(cvgnpz[7,7]);
stdzc=sqrt(cvgnpz[9,9]);
stdzx=sqrt(cvgnpz[10,10]);
stdgnp=sqrt(cvgnpz[1,1]);
sdti=sqrt(cvgnpz[21,21]);
stdl=sqrt(cvgnpz[22,22]);
stdpy=sqrt(cvgnpz[6,6]);
stdpx=sqrt(cvgnpz[5,5]);
                @relative std deviations@

rstdx=stdx/stdgnp;
rstdc=stdcs/stdgnp;
rstdlc=stdlc/stdgnp;
rstdkx=stdkx/stdgnp;
rstdzc=stdzc/stdgnp;
rstdy=stdy/stdgnp;
rstdi=sdti/stdgnp;
rstdl=stdl/stdgnp;


            @Printing results@

print"
--------------------------- Standard Deviations ------------------------------";
print "        Output          C               X                  I                L               Y";
format/rd 12,4;
print "";
print 1~rstdc~rstdx~rstdi~rstdl~rstdy;
print "           (-)            (.73)            ( -)              (3.20)         (1.16)           (-)
";
print"
------------------------ Correlations with Output ----------------------------";
print "      Output            C               X                 I                  L               Y"
;
format 12,4;
print "";
print 1~cv[1,2]~cv[1,3]~cv[1,21]~cv[1,22]~cv[1,4];
print "         (-)              (.82)             (-)             (.90)            ( .86 )             (-)
";
print"
--------------------------- AR(1) Coefficients -------------------------------";
print "Output       C             X              I              L                Y";

format 3,2;
print arcv[1,1]~arcv[2,2]~arcv[3,3]~arcv[21,21]~arcv[22,22]~arcv[4,4];

print "(.90)         (.84)          (-)            (.76)       ( .90 )            (-)";




endif;

						@SIMULATIONS@

lcs=lc/(lc+lx+ly);
lys=ly/(lc+lx+ly);
lxs=lc/(lc+lx+ly);



if rows(inx)==1;
yinp=0~0~0~0~0~0~0~0~0;
t=1;
bigt=500;bigt=50;

y=0|0|0|0|0|0;
y=yt|0|0|0|0|0|0;




do while t<bigt;sx=rndn(1,1)*swx;sy=rndn(1,1)*swy;ep1=rndn(1,1)*ewc;
ep2=rndn(1,1)*ewx;ep3=rndn(1,1)*ewy;

@next five lines for impulse response@

If t==1;
sx=0.00039;sy=0.00039;ep1=0.00039;ep2=0.00039;ep3=0.00039;
else;
sx=0.000;sy=0.000;ep1=0.000;ep2=0.000;ep3=0.000;
endif;

@Note: yt=khx,khy,phx,zhc,zhx,zhy@
yt=nj*yt+qqdd*(0|0|sx|ep1|ep2|ep3);

ny=yt[3,.]|(vdy*yt)|yt[3,.]|(vdy*yt)|yt[3,.]|(vdy*yt)|yt[1:2,.]|0;

inp=ss*ny;
nch=(al0+a0)*inp[7]+(al1+a1)*inp[4]+(al2+a2)*inp[1]+yt[4,1];
nxh=(be0+b0)*inp[8]+(be1+b1)*inp[5]+(be2+b2)*inp[2]+yt[5,1];
nyh=(ga0+c0)*inp[9]+(ga1+c1)*inp[6]+(ga2+c2)*inp[3]+yt[6,1];
yinp=yinp|inp';
@NOTE:y=khx,khy,phx,phy,zhc.zhx,zhy,ch,phx+xh,phy+yh,gnp,l@

y=y~(yt[1:3,.]|(vdy*yt)|yt[4:6,.]|nch|(yt[3,.]+nxh)*sixs|((vdy*yt)+nyh)*siys|
xgnp*(yt[3,.]+nxh)+((vdy*yt)+nyh)*ygnp+cgnp*nch|lcs*inp[7]+lxs*inp[8]+lys*inp[9]);
t=t+1;
endo;

else;
t=1;
yt=0|0|0|0|0;
y=yt|0|0|0|0|0|0|0;


bigt=500;bigt=	150;
	

do while t<bigt;sx=rndn(1,1)*swx;sy=rndn(1,1)*swy;ep1=rndn(1,1)*ewc;
ep2=rndn(1,1)*ewx;ep3=rndn(1,1)*ewy;

				@next five lines for impulse response@


If t==1;
sx=0.00039;sy=0.00039;ep1=0.00039;ep2=0.00039;ep3=0.00039;
else;
sx=0.000;sy=0.000;ep1=0.000;ep2=0.000;ep3=0.000;
endif;

@Note: yt=khx,khy,zhc,zhx,zhy@
yt=nj*yt+qqdd*(0|0|ep1|ep2|ep3);

ny=(vvr*yt)|(vvr*yt)|(vvr*yt)|yt[1:2,.]|0;
ppp=vvr*yt;
inp=ss*ny;
nch=(al0+a0)*inp[7]+(al1+a1)*inp[4]+(al2+a2)*inp[1]+yt[3,1];
nxh=(be0+b0)*inp[8]+(be1+b1)*inp[5]+(be2+b2)*inp[2]+yt[4,1];
nyh=(ga0+c0)*inp[9]+(ga1+c1)*inp[6]+(ga2+c2)*inp[3]+yt[5,1];

@NOTE:y=khx,khy,phx,phy,zhc.zhx,zhy,ch,phx+xh,phy+yh,gnp,l@

y=y~(yt[1:2,.]|(vvr*yt)|yt[3:5,.]|nch|(nxh+ppp[1])*sixs
|(ppp[2]+nyh)*siys|
(nxh+ppp[1])*xgnp+(ppp[2]+nyh)*ygnp+nch*cgnp|lcs*inp[7]+lxs*inp[8]+lys*inp[9]);

t=t+1;

endo;
endif;
y=y';
				@GRAPHING@

y=real(y);   @Even though the numbers are real, gauss thinks they are complex. This messes up 
		      the maxc and minc commands, which use the modulus, and yields absolute values. 
                      But maxc and minc are used by xy in sizing the window and plotting. This messes up 
	              the graph. Using the "real" command fixes the bug.@

tt=seqa(1,1,bigt);
zz=zeros(bigt,1);

graphset;					@ reset graph parameters        @

_pmcolor=15;_pltype={1 3 5 6};_plsel={1 2 4 5 6};_ptek="kaz10.tkf";
                                   
_plegctl=1;                               @ set loc. and size of legend   @
_plegstr= "y\000i\000c";                  @ legends @
xlabel("time");                             @ x-axis label                  @
titstr="Impulse Response for GNP, C, and I";         @  graph title                   @
@titstr="Simulated Data";            graph title                   @
title(titstr);                              @ put title                     @
_pdate=0;                                   @ no date                       @
@_pltype = {6 1 3};                           line types                  @


@graphprt("  -pr=3 -c2 -cf =c:\\data\\impres.pic  -CL=L -cs=h ");@



xy(tt,y[.,11]~(y[.,9]+y[.,10])~y[.,8]);

@stdc(y[.,11]~(y[.,9]+y[.,10])~y[.,8]);@
end;