这个是程序(分为了两部分):
/*
@>>>>>>>>>>>>>>>>>>>>Time varying parameter model >>>>>>>>>>>@
@>>>>>>>>>>>>>>>>>>>>with Markov Varance >>>>>>>>>>>>>>>>>>>@
@>>>>>>>>>>>>>>>>>>>>Using Optmum Proc >>>>>>>>>>>>>>>>>>>>>>@
@ May 9, 1990 version @
cjkim@korea.ac.kr
*/
new;
library optmum,pgraph;
load data[120,6]=tvpmrkf.prn;
/* column 1: m1===growth rate of quarterly average M1
2: dint=change in the lagged interest rate (3-month T-bill)
3: inf==lagged inflation
4: surp==full employment budget surplus (detrended)
5: m1lag==lag of m1
6: quarter index
1959.3--1987.4, ........................ */
qtr_idx=data[.,6];
T=rows(data);
yy=data[.,1];
xt=ones(t,1)~data[.,2:5];
output file=tvpmrkf.out reset;
output off;
@================= Initialize Global Variables============@
START=11; @1962.1--1989.2 @
PRMTR_NO=9; @ >>Number of parameters to be estimated. >>@
TVP_NO=5; @ >>Number of Time-varying Coefficients >>@
pr_1=sqrt(0.2); pr_2=sqrt(0.2); pr_3=sqrt(0.2);
pr_4=sqrt(0.2);pr_5=sqrt(0.2);pr_6=sqrt(0.2);
pr_7=sqrt(1.050); pr_8=0.9;pr_9=0.9;
@initial values of parameters@
PRMTR_IN=pr_1|pr_2|pr_3|pr_4|pr_5|pr_6|pr_7|pr_8|pr_9;
prmtr_in=0.36072|0.07080|0.079|0.0791|0.070|0.072|0.89439|
0.93266|0.95990;
output on; "initial values of prmtr is(vre,intcpt,dint,inf,surp,mlag,w,p,q";
prmtr_in';
output off;
PRIOR_CF=zeros(TVP_NO,1); @Initial values of reg coefs@
prior_vr=eye(tvp_no)*100;
@ Maximum Likelihood Estimation @
@==================================================@
{xout,fout,gout,hout}=optmum(&lik_fcn,tran_inv(PRMTR_IN));
@ prmtr estimates, -log lik value, Grandient, Hessian@
"Calculating Hessian..... Please be patient!!!!";
hout0=hessp(&lik_fcn,xout);
hout=inv(hessp(&lik_fcn,xout));
PRM_FNL=TRANS(xout); @ Estimated coefficients, constrained@
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@
output on;
"==FINAL OUTPUT========================================================";
"likelihood value is ";; fout;
"Estimated parameters are:";
prm_fnl';
"var-cov matrix is:"; Hsn_fnl;
"Standard errors of parameters are:"; sd_fnl;
"===============================================================";
OUTPUT OFF;
F_SS=SAVE_OUT(xout);
OUTPUT FILE=TVPMRKF.DTA RESET;
F_SS; OUTPUT OFF;
xy(seqa(1,1,rows(f_ss)),f_ss[.,1]);
xy(seqa(1,1,rows(f_ss)),f_ss[.,2:4]);
END;
@========================================================================@
@========================================================================@
PROC LIK_FCN(PRMTR1);
local prmtr, ppr,qpr,aaa,prob__0,prob__1,p_cf_0,p_cf_1, p_vr_0,
p_vr_1, qq, lik, j_iter, h, pr__0_l,pr__1_l, pst_cf0, pst_cf1,
pst_vr0, pst_vr1, f_cast0, f_cast1, ss00, ss01, ss10, ss11,
pr_vl00, pr_vl01, pr_vl10, pr_vl11, pr_val,k_gn00, k_gn01,
k_gn10, k_gn11, p_cf00, p_cf01, p_cf10, p_cf11, p_vr00, p_vr01,
p_vr10, p_vr11, pro_00, pro_01, pro_10, pro_11, likv;
EXTERNAL PROC TRANS, V_PROB;
PRMTR=TRANS(PRMTR1);
PPR=PRMTR[8,1]; @Pr[St=1/St-1=1]@
QPR=PRMTR[9,1]; @Pr[St=0/St-1=0]@
AAA=EYE(TVP_NO);
@>>>>>>>>>>>>>>>>>>>>>>>>> INITIAL PROB. Pr[S0/Y0] @
PROB__1=(1-QPR)/(2-PPR-QPR); @Pr[St-1=1/Yt-1], STEADY STATE PROB.@
PROB__0=1-PROB__1 ; @Pr[ST-1=0/Yt-1], STEADY STATE PROB @
P_CF_0=PRIOR_CF; P_CF_1=PRIOR_CF;
P_VR_0=PRIOR_VR; P_VR_1=PRIOR_VR; @ Initial values@
QQ=(PRMTR[2,1]^2~0~0~0~0)|
(0~PRMTR[3,1]^2~0~0~0)|
(0~0~PRMTR[4,1]^2~0~0)|
(0~0~0~PRMTR[5,1]^2~0)|
(0~0~0~0~PRMTR[6,1]^2);
LIKV=0.0;
J_ITER = 1;
DO UNTIL J_ITER>T;
H=xt[J_ITER,.];
PR__0_L= QPR*PROB__0 + (1-PPR)*PROB__1; @Pr[St=0/Yt-1]@
PR__1_L= (1-QPR)*PROB__0 + PPR*PROB__1; @Pr[St=1/Yt-1]@
@===================PREDICTION=============================@
PST_CF0 = AAA * P_CF_0; @ Prediction of coef when St-1=0 @
PST_CF1=AAA*P_CF_1;
PST_VR0 = AAA * P_VR_0 * AAA' +QQ;
PST_VR1 = AAA * P_VR_1 * AAA' +QQ;
F_CAST0 = yy[J_ITER,.] - H * PST_CF0; @ Forecast error when St-1=0 @
F_CAST1 = yy[J_ITER,.] - H * PST_CF1; @ Forecast error when St-1=0 @
SS00= H * PST_VR0 * H' + PRMTR[1,1]^2; @ St-1=0,St=0 @
SS01= H* PST_VR0 * H' + PRMTR[7,1]^2; @ St-1=0,St=1 @
SS10= H * PST_VR1 * H' + PRMTR[1,1]^2; @ St-1=1,St=0 @
SS11= H * PST_VR1 * H' + PRMTR[7,1]^2; @ St-1=1,St=1 @
@<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<@
PR_VL00=V_PROB(F_CAST0,SS00)* QPR*PROB__0; @Pr[St,Yt/Yt-1]@
PR_VL01=V_PROB(F_CAST0,SS01)*(1-QPR)*PROB__0;
PR_VL10=V_PROB(F_CAST1,SS10)*(1-PPR)*PROB__1;
PR_VL11=V_PROB(F_CAST1,SS11)* PPR*PROB__1;
PR_VAL=PR_VL00+PR_VL01+PR_VL10+PR_VL11; @Pr[Yt/Yt-1]@
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
LIK=-1*LN(PR_VAL);
@=======================UPDATING===============================@
K_GN00 =PST_VR0 * H' / SS00;
K_GN01 =PST_VR0 * H' / SS01;
K_GN10 =PST_VR1 * H' / SS10;
K_GN11 =PST_VR1 * H' / SS11;
P_CF00 = PST_CF0 + K_GN00 * F_CAST0;
P_CF01 = PST_CF0 + K_GN01 * F_CAST0;
P_CF10 = PST_CF1 + K_GN10 * F_CAST1;
P_CF11 = PST_CF1 + K_GN11 * F_CAST1;
P_VR00 = (EYE(TVP_NO) - K_GN00 * H ) * PST_VR0;
P_VR01 = (EYE(TVP_NO) - K_GN01 * H ) * PST_VR0;
P_VR10 = (EYE(TVP_NO) - K_GN10 * H ) * PST_VR1;
P_VR11 = (EYE(TVP_NO) - K_GN11 * H ) * PST_VR1;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PRO_00=PR_VL00/PR_VAL; @Pr[St,St-1/Yt]@
PRO_01=PR_VL01/PR_VAL;
PRO_10=PR_VL10/PR_VAL;
PRO_11=PR_VL11/PR_VAL;
PROB__0=PRO_00+PRO_10; @Pr[St=0/Yt]@
PROB__1=PRO_01+PRO_11; @Pr[St=1/Yt]@
@<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<@
P_CF_0=(PRO_00*P_CF00 + PRO_10*P_CF10)/PROB__0;
P_CF_1=(PRO_01*P_CF01 + PRO_11*P_CF11)/PROB__1;
P_VR_0=(PRO_00*(P_VR00+(P_CF_0-P_CF00)*(P_CF_0-P_CF00)')+
PRO_10*(P_VR10+(P_CF_0-P_CF10)*(P_CF_0-P_CF10)'))/PROB__0;
P_VR_1=(PRO_01*(P_VR01+(P_CF_1-P_CF01)*(P_CF_1-P_CF01)')+
PRO_11*(P_VR11+(P_CF_1-P_CF11)*(P_CF_1-P_CF11)'))/PROB__1;
IF J_ITER < START; goto skip; endif;
LIKV = LIKV+LIK;
skip: J_ITER = J_ITER+1;
ENDO;
RETP(LIKV);
ENDP;