/* nonling.src: Collection of procedures to execute likelihood evaluation using the particle filter.
		The likelihood function is associated with the optimal growth model, outlined in Section 10.2.2.2
*/

/* The collection of procedures contained in ara.prg facilitate the construction of the policy function for determining
	consumption. Inputs provided by these procedures are denoted by lines ending with //*** 
*/

/* Begin with a collection of globals */

nparts = 20000;		// size of particle swarms
lnlkadj = 400;		// adjustement to the likelihood function that prevents underflows
big0 = zeros(nparts,1);	// vector that will be used to execute resampling 
big1 = ones(nparts,1);	// vector that will be used to execute resampling 

//common random numbers
rndseed 12345678;
//use for drawing of initial state of a and k in partproc
crn_rndn0=rndn(2,nparts);
//use for epsilons in partproc
crn_rndn1=rndn(nobs,nparts);
//use for resampling in partproc
crn_rndu0=rndu(nobs,nparts);


proc (3) = lkeval(datt,veet,s_1,opt_gam_from_polcfn,wmat,lnnormcons);

    /* inputs: 
		datt: 	time-t data (y,i)
		veet: 	time-t structural shocks (eps)
		s_1:	time-(t-1) state (a,k) 
		opt_gam_from_polcfn:	policy function parameters //***
		wmat:	inverse of vcv matrix of errors 
		lnnormcons: normalizing constant for likelihood function: ln( ((2*pi)^(-1.0))*( (det(wmat))^0.5 ) )
    */

    /* output: at, kt, likelihood value */

    local z_1,k_1,eps,y,i,ztilde,ktilde,c_1,y_1,i_1,k,z,ypred,cpred,ipred,uy,ui,u,
            lnkern,lnlkvaladj,lkvaladj,;

    z_1 = s_1[1];
    k_1 = s_1[2];
    eps = veet;
    y = datt[1];
    i = datt[2];

    /* convert yesterday's a and k into today's a and k. 
	the 'if' statements check that states lie between lower and upper bounds */

    if k_1<kl;
        k_1=kl;
    elseif k_1>kh;
        k_1=kh;
    endif;
    if z_1<zl;
        z_1=zl;
    elseif z_1>zh;
        z_1=zh;
    endif;

    ztilde = z_1/omegaz;		// normalizes z to lie between [-1 1]
    ktilde = (k_1-kss)/omegak;		// normalizes k to lie between [-1 1]


    c_1 = ceeofzk(ztilde,ktilde,opt_gam_from_polcfn);	//***
    y_1 = exp(z_1)*(k_1^alp);		// alp = alpha
    i_1 = y_1-c_1;
    k = i_1 + (1 - del)*k_1;		// del = delta
    z = rho*z_1 + eps;  

    /* construct predictions for todays y, c, i */

    if k<kl;
        k=kl;
    elseif k>kh;
        k=kh;
    endif;
    if z<zl;
        z=zl;
    elseif z>zh;
        z=zh;
    endif;

    ztilde = z/omegaz;
    ktilde = (k-kss)/omegak;

    ypred = exp(z)*(k^alp);
    cpred= ceeofzk(ztilde,ktilde,opt_gam_from_polcfn);	//***
    ipred = ypred-cpred;

    /* construct implied errors */

    uy = y - ypred;
    ui = i - ipred;
    u = uy|ui;
    lnkern = -0.5*(u'wmat*u);
    lnlkvaladj = lnnormcons+lnkern+lnlkadj; @ an adjustement term (lnlkadj) is added 
                                                to prevent underflow; eliminated below @ 
    lkvaladj = exp(lnlkvaladj); 

    retp(z,k,lkvaladj);

endp;

proc partproc(opt_gam_from_polcfn,wmat,lnnormcons,cvcvs0tilde);

	local rsnewz,rsnewk,rslk,lts,lnlts,jjj,partz,partk,parteps,partuy,partui,newz,newk,lk,iii,
		s0tildedraw,probs,csprobs,udraw,inddiff,drawindvec,drawind,lnlk


	;

   /* inputs: 
		opt_gam_from_polcfn:	policy function parameters //***
		wmat:			inverse of vcv matrix of errors 
		lnnormcons: normalizing constant for likelihood function: ln( ((2*pi)^(-1.0))*( (det(wmat))^0.5 ) )
		cvcvs0tilde:		vcv matrix of s0
*/

    /* output: ln likelihood value */


rsnewz = zeros(1,nparts);    /* will contain resampled updated (over time) as and ks: a' and k' */
rsnewk = zeros(1,nparts);
rslk = zeros(1,nparts);     /* will contain resampled likelihood values */

lts = zeros(nobs,1);         /* will contain date-t (averaged) likelihood values */
lnlts = zeros(nobs,1);      /* will contain date-t log likelihood values */

jjj = 1; do while jjj<=nobs; /* index over dates */

partz = zeros(1,nparts);        /* will contain initial sample of as and ks, epsilons, uy's, uc's, ui's */
partk = zeros(1,nparts);
parteps = zeros(1,nparts);
partuy = zeros(1,nparts);
partui = zeros(1,nparts);
newz = zeros(1,nparts);         /* will contain updated as and ks pzssed to the next round */
newk = zeros(1,nparts);
lk = zeros(1,nparts);

iii = 1; do while iii<=nparts;  /* index over particles */

    /* draw initial state */

    if jjj == 1;

        s0tildedraw = cvcvs0tilde*crn_rndn0[.,iii];
        partz[1,iii] = s0tildedraw[2];
        partk[1,iii] = kss*exp(s0tildedraw[1]);

    	parteps[1,iii] = crn_rndn1[jjj,iii]*sige; 	// sige is the standard deviation of epsilon
                        
    else;
        partz[1,iii] = rsnewz[1,iii];
        partk[1,iii] = rsnewk[1,iii];
	    parteps[1,iii] =    crn_rndn1[jjj,iii]*sige; 
                            
    endif;
    { newz[1,iii],newk[1,iii],lk[1,iii] } = 
		lkeval(rawdat[jjj,.]',parteps[1,iii],partz[1,iii]|partk[1,iii],opt_gam_from_polcfn,wmat,lnnormcons);
			// rawdat contains the data

 /* newa and newk are updated values: a' and k' */

iii = iii+1; endo;
lk = normal(lk);		//prevents underflow: see the proc normal below
probs = lk./(sumc(lk'));
csprobs = cumsumc(probs');

lts[jjj] = meanc(lk');
lnlts[jjj] = ln(lts[jjj])-lnlkadj;  @ the term -lnlkadj eliminates the adjustement term @

/* resampling */

iii = 1; do while iii<=nparts;  /* index over particles */

    udraw = crn_rndu0[jjj,iii].*big1;
    inddiff = (csprobs-udraw) .>= big0;
    drawindvec = indexcat(inddiff,1);
    drawind = drawindvec[1];  

    /* draw new state */

    rsnewz[1,iii] = newz[1,drawind];
    rsnewk[1,iii] = newk[1,drawind];
    rslk[1,iii] = lk[1,drawind];

iii = iii+1; endo;

jjj = jjj+1; endo;

lnlk = sumc(lnlts);

retp(lnlk);

endp;

proc normal(x);
    retp(x.*(abs(x) .> 4.19e-307));
endp;
