/* 
** A Time-Varying Parameter Model of U.S. Monetary Growth Function 
** (Kim and Nelson, JBES 1989)
** DM = b0 + b1 DR1 + b2 DP1 + b3 SURP1 +b4 DM1 + u
** b(t) = b(t-1) + v(t), where b={b0,b1,b2,b3,b4}
** u(t) ~ nii(0,R), R is 1x1 (r1)
** v(t) ~ nii(0,Q), Q is 5x5 assuming zeros off-diagnols
** (q0, q1, q2, q3, q4 corresponds to b's)
*/
use gpe2;
n=107;
load data[n,6]=c:\course06\XMU\tvp.txt; @ Kim and Nelson (1989) TVP.DAT @
/* Data description: quartely 1959.3-1985.4
   1: Quarter Index
   2: dm    = growth rate of quarterly average M1
   3: dr1   = change in the lagged interest rate (3-month T-bill)
   4: dp1   = lagged inflation
   5: surpl = lagged full employment budget surplus
   6: dm1   = lag of m1 (1959.2--1985.3)
** Data used for estimation: 1962.1 - 1985.4 (Obs.11-106)
   Number of time-varying parameters = 5
   Number of parameters to be estimated = 6
   Initial values (Kim and Nelson)
   _b = {0.5,0.1,0.1,0.1,0.1,0.1}
*/
qtr=data[2:n,1];
data=data[2:n,2:6];

call reset;
_nlopt=1;
_method=5;
_conv=1;
_iter=100;
b0={0.5,0.1,0.1,0.1,0.1,0.1};

call optimize(&kf_llf,data,b0);

b=kf_state(data,__b);

print b;

end;

/* 
** user defined kalman filter model for each observation
** (from z and b, compute paramters for state-space model)
** because some of state-space parameters may be data dependent, 
** that is, they are different for different observation. 
** kalman filter is defined for each observation, one at a time. 
** inputs: z = observed data matrix, b = parameter vector, i = obs 
** outputs: y (may be more than 1) = endogenous variable(s),
** h, a, c, f, r, q = state-space parameters
*/
proc (7) = kf_model(z,b,i);
    local y,x,h,a,c,f,r,q;
    y=z[i,1];     @ endogenous variable(s) @
    x=1~z[i,2:5]; @ predetermined variables incl. constant @

    a=0;
    h=x;
    r=b[1]^2;
    c=zeros(5,1);
    f=eye(5);
    q=diagrv(eye(5),b[2:6]^2);

    retp(y,h,a,c,f,r,q);
endp;

/* 
** initial values of state variables and var-cov matrix
** in general, use the unconditional mean and covariance
** for the case of stationary state variables.
** if the state variables are nonstationary, the initial
** values can be set arbitrarily, e.g., zeros vector.
** then covariance matrix must have large values to the 
** diagonal elements.
*/
proc (2) = kf_initial(z,b);
    local y,h,a,c,f,r,q,k,b0,p0;

    {y,h,a,c,f,r,q}=kf_model(z,b,1);
    k=rows(c); @ rows(c)=rows(f)=cols(f)=cols(h) @
    /*
    b0=inv(eye(k)-f)*c;
    p0=reshape(inv(eye(k*k)-f.*.f)*vec(q),k,k)';
    */
    @ wild guess of b0 should come with large p0 @
    b0=zeros(k,1);
    p0=eye(k)*50;

    retp(b0,p0);
endp;

/*
** KALMAN1.GPE
** Kalman Filter
** (1) Y(t)=H(t)b(t)+a(t)+u(t) ~ nii(0,R(t))
** (2) b(t)=F(t)b(t-1)+c(t)+v(t) ~ nii(0,Q(t))
** Y: NxM data matrix, Y(t): Mx1
** H(t): MxK paramter matrix 
** a(t): Mx1 parameter vector 
** b(t): Kx1 parameter vector (unobsered state variables)
** F(t): KxK parametr matrix (first-order transition matrix)
** c(t): Kx1 parameter vector 
** R(t): MxM paramter matrix (var-cov matrix of (1))
** Q(t): KxK parameter matrix (var-cov matrix of (2))
** a,H,R,c,F,Q may be constant or data dependent on r exogenous or 
** predetermined data variables X(t) (X: Nxr data matrix) 
**
** To use this module, the user needs to define the model (in the
** procedure kf_model) and initial values for state variables 
** (in the procedure kf_initial). log-likeligood function kf_llf
** and state variable inference kf_state are compued based on 
** kf_model and kf_intital.
*/

/*
** log-likelihood function based on 
** kf_model and kf_initial
*/
proc kf_llf(z,b);
    local n,y,h,a,c,f,r,q,b0,p0,b1,p1;
    n=rows(z);    
    {b0,p0}=kf_initial(z,b);
    
    local i,e,v,g,ll;
    ll=0;
    i=1;
    do until i>n;
        {y,h,a,c,f,r,q}=kf_model(z,b,i);

        /* prediction */
        b1=c+f*b0;
        p1=f*p0*f'+q;

        @ error and variance @
        e=y'-h*b1-a; 
        v=h*p1*h'+r;

        @ log-likelihood @
        if i>10;
        ll=ll-0.5*cols(y)*ln(2*pi)-0.5*ln(det(v))-0.5*e'*invpd(v)*e;
        endif;

        /* updateting */
        @ kalman gain @
        g=p1*h'*invpd(v);

        b0=b1+g*e;
        p0=p1-g*h*p1;

        i=i+1;
    endo;

    retp(ll);

endp;

/*
** state variable inference: basic filter and smoothing
** given the estimated parameters, predicted the state variables
** based on kf_model and kf_initial
*/
proc kf_state(z,b);
    local n,y,h,a,c,f,r,q,b0,p0,b1,p1;
    n=rows(z);

    {b0,p0}=kf_initial(z,b);

    /* save state variables (and covariances) */
    local k,beta,beta1,betas,omega,omega1,omegas;
    k=rows(b0);
    beta=miss(zeros(n,k),0);
    omega=miss(zeros(n,k*k),0);
    beta1=beta;
    omega1=omega;
    @ beta (omega): basic filtered paramters (covariances) @
    @ betas (omegas): smooth filtered parameters (covariances) @

    local i,e,v,g,ll;
    ll=0;
    i=1;
    do until i>n;
        {y,h,a,c,f,r,q}=kf_model(z,b,i);

        /* prediction */
        b1=c+f*b0;
        p1=f*p0*f'+q;

        @ error and variance @
        e=y'-h*b1-a; 
        v=h*p1*h'+r;

        /* updateting */
        @ kalman gain @
        g=p1*h'*invpd(v);

        b0=b1+g*e;
        p0=p1-g*h*p1;

        @ keep predicted state variables and covariances @
        if i>10; 
            beta[i,.]=b0';
            omega[i,.]=vec(p0)';
            beta1[i,.]=b1';
            omega1[i,.]=vec(p1)'; 
        endif;

        i=i+1;
    endo;
    
    @ smoothing @
    betas=beta;
    omegas=omega;
    i=n-1;
    do until i<10;
        {y,h,a,c,f,r,q}=kf_model(z,b,i+1);
        p0=reshape(omega[i,.],k,k)';
        p1=reshape(omega1[i+1,.],k,k)';
        g=p0*f'*invpd(p1);
        betas[i,.]=beta[i,.]+(g*(betas[i+1,.]'-c-f*beta[i,.]'))';
        omegas[i,.]=vec(p0+g*(reshape(omegas[i+1,.],k,k)'-p1)*g')';
        i=i-1;
    endo;
   /*
    print beta;
    print omega;

    print betas;
    print omegas;
    */

    /* plots of beta or betas (and s.e.) */ 
    pqgwin many;
    library pgraph;
    graphset;
    i=1;
    do until i>k;
        b0=betas[.,i];
        p0=sqrt(omegas[.,(i-1)*k+i]); @ s.e. of betas[.,i] @
        call xy(seqa(1,1,n),b0~(b0+2*p0)~(b0-2*p0));
        i=i+1;
    endo;
    
    retp(betas);
endp;