/*  stogrow.src:  Procedures designed to solve the stochastic optimal growth model */

@  
Ultimate goal: convert system of nonlinear model equations into a linear system
of the form

(1)	x(t) = Fx(t-1) + Ge = Fx(t-1) + v(t), 

where Evv' = Q,
x(t) is an nx1 vector of model variables,
and e is a jx1 vector of exogenous disturbances.

Observable variables are presumed to be contained in the jx1 vector
y(t) (note there are as many observable variables as exogenous disturbances).
The mapping between model variables and observables is given by

(2)	y(t) = H'x(t), where H is an nxj observer matrix.

The matricies Y (Txj matrix of observed values of y(t), t = 1,...,T), F, H, and Q are
passed to the Kalman filter procedure for likelihood evaluation.

Prior to obtaining (1), we take the intermediate step of obtaining

(3)	G0x(t) + G1x(t-1) = GEe(t),

where Gj is nxn, j = 0, 1, with the ith row corresponding to the ith log-linearized 
model equation.  The mx1 vector e(t) contains the variables v(t), as well as expectational
errors associated with the euler equations of the model.  The matrix GE is nxm.  It maps
the appropriate elements of e(t) into the appropriate linearized model equations.

Model parameters are contained in the kx1 vector p.  Ordering: 

p = [alp|bet|del|phi|rho|sige]'

Logs of variables are contained in the nx1 vector x.  Ordering:
x = [y c i k a ]';

Logged steady state values of x are in xbar.

The procedures:

--- steady(p):  converts parameter values into steady state values of x.

---  loglin(xbar):  evaluates logs of the k model equations at xbar.  When serving as input
to gradp, computes log-linear approximations of each equation.  The format for
calling this proc is as follows:

G = gradp(&loglin,xbar);

Output is the nxn matrix G, whose ith row contains partial derivatives of the log of
the ith model equation with respect to the log ratio of each variable to its steady
state value.  The purpose of this is to convert the nonlinear model into a linear system
of the form (2).  When the variable tme is set to j, gradp(&loglin,xbar) returns Gj, j = 0, 1.

---  qproc(G0,GE): creates the matrix Q from G0 and GE.  returns q and a version of 
				g called gno under which the shocks are uncorrelated.  this can be used
				to construct impulse responses to uncorrelated shocks in endowsts.src

---  hproc(p):  converts model parameters into the observer matrix

--- modelsol(p): reads in p, generates as output the system matricies F, GeNO, Q and H.  Note
				that GeNO can be used to generate impulse responses to UNCORRELATED
				shocks

NOTE: modelsol is NOT model-specific!!
@


proc steady(p);
    /* calculates the steady state of the model */
    local alp,bet,del,kss,yss,iss,css;
    alp = p[1];
    bet = p[2];
    del = p[3];

    kss = (alp/(1/bet-1+del))^(1/(1-alp));
    yss = kss^alp;
    iss = del*kss;
    css = yss - iss;

    retp(yss|css|iss|kss|1);
endp;

proc loglin(x);
local y, ylag, c, clag, i, ilag, k, klag, a, alag, alp, bet, del, rh, sige, ph, cfac, lfac, z;

	alp = p[1];
	bet = p[2]; 
	del = p[3]; 
	rh  = p[5]; 
	sige= p[6]; 
	ph = p[4];

	@  
		Note: both current and lagged variables are set to logged ss values. When tme == 0,
		derivatives are taken with respect to current variables; when tme == 1, derivatives
		are taken with respect to lagged variables.
	@

	if tme ==0;                                                          @ Build G0 @
		y = x[1]; c = x[2]; i = x[3]; k = x[4]; a = x[5];
		ylag = xbar[1]; clag = xbar[2]; ilag = xbar[3];  
		klag = xbar[4]; alag = xbar[5];
	else;                                                                     @ Build G1 @
		y = xbar[1]; c = xbar[2]; i = xbar[3];  
		k = xbar[4]; a = xbar[5];
		ylag = x[1]; clag = x[2]; ilag = x[3];  
		klag = x[4]; alag = x[5];
	endif;


	z = zeros(nvars,1);

	z[1] = -ph*clag - ln(bet) + ph*c - ln( alp*exp(a)*exp( (alp-1)*klag ) + (1-del) );
	z[2] = y - a - alp*klag;
	z[3] = y - ln( exp(c) + exp(i) );
	z[4] = k - ln( exp(i) + (1-del)*exp(klag) );
	z[5] = a - rh*alag;

	retp(z);

endp;

proc (2) = qproc(p,g0inv);

	/* return gno for uncorrelated shocks; g for correlated shocks */
	
	local sige, sig, cov, q, gee, geeno, g, gno;
	
	sige = p[6];
	sig = zeros(nshocks,nshocks);
	gee = zeros(nvars,nvars);
	geeno = zeros(nvars,nvars);
	sig[1,1] = sige^2; 
	cov = chol(sig)';
	gee[5,5] = cov;
	geeno[5,5] = diagrv(eye(nshocks),sqrt(diag(sig)));
	g = g0inv*gee;
	q = g*g';
	gno = g0inv*geeno;

	retp(g,q);

endp;

proc hproc(p);

	local hmat;
	hmat = zeros(nvars,nshocks);
	hmat[1,1] = 1;

	retp(hmat);

endp;

proc (6) = modelsol(p);                      @ Takes as input a vector of parameters. Generates
                                                                       as output the system matricies F, Q and H @

local xbar, g0, g1, ng0, ng1, gmat, va, vi, cmat, vei, i, invcmat, invng0, fmat, hmat,
qmat, gnomat;

xbar = ln(steady(p));  
tme = 0;
g0 = gradp(&loglin,xbar);
tme = 1;
g1 = gradp(&loglin,xbar); 
g1 = -1*g1;

if rank(g0) < nvars;
	retp("Noninvertible g0",0,0,0);
	goto endd;
endif;

gmat = inv(g0)*g1;
{va,vi,cmat,vei} = eigrg2(gmat');        @  va: eigenvalues.  cmat: associated eigenvectors @

i = indexcat(abs(va),1|1000);    @  Returns row indexes indicating the location of eigs
                                                                          greater than 1 (and <= 1000)  @

if rows(i) /= neulers or ismiss(i);       @  Checks whether number of explosive
                                                                                     eigs equals number of Euler Eqns  @
	retp("wrong number of eigenvalues",rows(i),0,0);
	goto endd;
endif;

invcmat = cmat';
ng0 = g0;
ng1 = g1;
ng0[eulers,.] = invcmat[i,.];
ng1[eulers,.] = zeros(neulers,nvars);
invng0 = inv(ng0);

if ismiss(invng0);
	retp("non-invertible ng0",0,0,0);
	goto endd;
endif;

fmat = invng0*ng1;
{gnomat,qmat} = qproc(p,invng0);
hmat = hproc(p);

retp(g0,g1,fmat,gnomat,qmat,hmat);

endd:
endp;




