我想对AR(2)的phi2, 也就是y_{t-2}的系数构造一个95%的likelihood ratio confidence interval,最大似然估计采用exact MLE,具体用kalman filter来做
model: (y_t - mu)=phi1*(y_{t-1} - mu) + phi2*(y{t-2} - mu) + epsilon_t, epsilon_t ~ N(0, var),所以,参数theta=( mu, var, phi1, phi2)
因为采用exact MLE, 所以必须限制| phi1 |<2,我想固定住phi1, 然后用 qnewton来求likelihood的最大值,并且求出相应的phi2的估计,如果这个最大值超出-341.368,那么phi2就在95%置信区间内,但是程序提示error G0003 : Indexing a matrix as a vector,这个应该是 Kalman Filter 里面用列向量 phi1 来定义矩阵F出了问题,但我不知道为什么会有问题。。。
我已经花了2个小时来修改程序,可是还是搞不定,希望板上有大侠能帮忙解答哈。。。感激不尽!
程序提示错误如下(刚开始没看到第一行的错误提示,汗。。。) :
H:\coursework\S1_cousework\Applied Ecometrix\hw3\Q2\ii\LR CI(58) :
error G0003 : Indexing a matrix as a vector
Currently active call: lik_fcn [59] H:\coursework\S1_cousework\Applied Ecometrix\hw3\Q2\ii\LR CI
Stack trace:
lik_fcn called from d:\src\qnewton.src, line 178
_bfgs called from d:\src\qnewton.src, line 98
qnewton called from H:\coursework\S1_cousework\Applied Ecometrix\hw3\Q2\ii\LR CI, line 27
程序如下 (data见附件):
new;
cls;
library pgraph;
format /m1 /rd 9,3;
load data[256,1]=gdp.txt;
yy=100*ln(data);
dy=yy[2:rows(yy)]-yy[1:rows(yy)-1];
T=rows(dy); @sample size@
n=1;
phi1=seqa(-2+0.00001,0.01,399);
s=rows(phi1);
phi2_mat={};
lik_mat={};
do until n > s;
@MAXIMUM LIKELIHOOD ESTIMATION@
prmtr_in=zeros(3,1); @starting values for parameters@
prmtr_in[1]=meanc(dy); @Use the sample mean as starting value for mu@
prmtr_in[2]=ln((stdc(dy))^2); @Use sample variance of dy as starting value for var_e@
alpha=0.05; @level at which to set Type I error@
{xout,fout,gout,cout}=qnewton(&lik_fcn,prmtr_in);
prm_fnl=trans(xout); @final constrained point estimates@
lik=-fout; @final likelihood value@
if lik ge -341.368;
phi2_mat=phi2_mat|prm_fnl[3];
lik_mat=lik_mat|lik;
endif;
n=n+1;
endo;
end;
@========================================================================@
@========================================================================@
proc LIK_FCN(PRMTR1); /* Kalman Filter*/
local prmtr,mu,var,phi1,phi2,theta1,theta2,f,q,h,x_ll,p_ll,
llikv,j,x_tl,p_tl,eta_tl,f_tl,x_tt,p_tt,vP_ll,A,z,R,y;
prmtr=trans(prmtr1); @constrain parameters@
@PARAMETERS, AR(2), phi1 fixed@
mu=prmtr[1]; @unconditional mean@
var=prmtr[2]; @standard deviation of error term@
phi2=prmtr[3]; @ar coefficient@
@SET UP STATE SPACE FORM@
F= phi1[n]~phi2|
1~0;
Q= var~0|
0~0;
H= 1~0;
z=ones(T,1); @Let "z" denote a vector of ones to pick up a constant term in the observation equation@
A= mu; @coefficient on the constant term@
R=0; @no error in the observation equation@
y=dy; @Let "y" denote the name of a generic observation variable or vector, even though it is real GDP growth in this case@
@INITIAL STATE ESTIMATES@
X_LL= zeros(rows(F),1); @unconditional expectation of state vector@
vP_LL = inv( (eye(rows(F)^2)-F.*.F) )*vec(Q);
P_LL= reshape(vP_LL,rows(F),rows(F));
llikv = 0;
j=1;
do until j>T;
@KALMAN FILTER@
X_TL = F * X_LL;
P_TL = F* P_LL * F' + Q;
ETA_TL = y[j] - H * X_TL - A*z[j];
F_TL = H * P_TL * H' + R; @lambda_t@
X_TT = X_TL + P_TL*H'*INV(F_TL) * ETA_TL;
P_TT = P_TL - P_TL * H'*INV(F_TL) * H * P_TL;
@LIKELIHOOD FUNCTION@
llikv=llikv - 0.5*(ln(2*PI*F_TL) + (ETA_TL^2)/F_TL);
@UPDATE for NEXT ITERATION@
X_LL=X_TT;
P_LL=P_TT;
j=j+1;
endo;
retp(-llikv); @return negative of likelihood because qnewton is a minimization routine@
endp;
@========================================================================@
@PROCEDURE TO CONSTRAIN PARAMETERS@
proc trans(c0);
local c1,aaa,ccc;
c1 = c0;
c1[2]=exp(c0[2]); @positive variance@
aaa=phi1[n]./2;
ccc=(1-abs(aaa))*c0[3]./(1+abs(c0[3]))+abs(aaa)-aaa^2;
c1[3]=-1* (aaa^2+ccc);
retp(c1);
endp;