📄 maxseek
字号:
@ This GAUSS file reads in data, sets options, and calls numerical
optimization routine for numerical estimation of SWARCH model @
output file=junk.out reset;
format /m1 /ros 16,8;
/* ======================================================================= */
@ Alter this section so as to read in your data and set parameters @
@ Crisp data @
capt = 1331; @ capt is the sample size @
load y[capt,1] = crisp.w;
@--------- Adjust any of the following to control specification desired --- @
nk = 5; @ nk is the first observation to be used in
estimation @
ns = 3; @ ns is the number of primitive states @
pphi = 1; @ pphi is the number of lags in autoregression
for y; @
izz = 1; @ izz = 1 means pij parameterized as
th(ij)^2/sum(th(ij)^2)
izz = 2 means pij parameterized as pij @
ipm =3; @ ipm specifies order in which transition probs
are parameterized
ipm = 1 implies p11 and p22 estimated
ipm = 2 implies pij for i=1,..,n j=1,..,n-1
ipm = 3 user input code @
irs = 2; @ irs denotes the number of probabilities that
are restricted @
karch = 2; @ karch is the order of the ARCH process @
larch = 1; @ larch = 0 means no leverage effect,
larch = 1 means leverage effect @
barch = 2; @ barch = 0 means estimate initial variance by MLE
barch = 2 means don't need initial variance @
tarch = 1; @ tarch = 0 means normal
tarch = 1 means t distribution @
@------ Input initial values for parameters ------------ @
/* The order in which variables are represented is as follows
constant term in regression
autoregressive terms in regression
initial variance parameter
constant term in ARCH equation
ARCH params in state 1
next elements: when izz =0 these are the transition probs
when izz =1 these are params v(i,j) such that
p(i,j) = v(i,j)^2 / sum j v(i,j)^2
elements are ordered as p11, p22 when ns =2
ordered as p(1,1),p(2,1),...,p(ns,1),
p(1,2),...,p(ns-1,ns) when ns > 2
next (ns - 1) elements: factor squared residuals are divided by
to get non-switching ARCH
leverage parameter
degrees of freedom parameter for t distribution */
proc startval; @procedure to set starting values @
local th;
let th[13,1] =
0.35080500 0.25003200 -0.56817500
0.027988000 -0.11706300
11.4016 0.051132 10.765493 0.1208310
4.3512660 13.146594 0.42415600 -5.2199370 ;
/* let th[11,1] =
0.34775 0.25589 0.55992 0.06036 0.16047 0.051751
0.99195 0.99907 5.7857 0.44348 3.7547 ; */
retp(th);
endp;
/*==========================================================================*/
@ In general no parts of this section should be changed @
ps = karch; @ ps is the number of lagged states that matter @
n = ns^(ps+1); @ n is the dimension of the state vector @
kc = 1; @ kc = 2 to print out ergodic probs and likelihood value
with each call to likelihood @
ks = 1; @ ks = 2 if smoothed probs are to be calculated @
captst = capt - nk +1; @ captst is the effective sample size @
resid = zeros(capt,1); @ regression residuals (filled in by procedures)@
skif = zeros(captst,n); @ skif is the matrix of filtered probs @
skis = zeros(captst,n); @ skis is the matrix of smoothed probs @
varfor = zeros(capt,1); @varfor is the vector of forecast variances @
outprob = zeros(capt,1); @outprob is the vector of outlier probabilities @
id = eye(ns); @ used in certain calculations below @
proc pattern1; @ This proc returns a (ps+1)*ns x n matrix. The ith
column contains a one in row j if st = j, contains a
one in row ns+j if st-1 = j, and so on @
local i1,ix,iq,na;
na = n/ns;
ix = eye(ns).*.ones(1,na);
i1 = 1;
do until i1 > ps;
na = na/ns;
iq = ones(1,ns^i1).*.(eye(ns).*.ones(1,na));
ix = iq|ix;
i1 = i1+1;
endo;
retp(ix);
endp;
/* ======================================================================= */
proc matpm(xth); @This proc defines the user's conventions for reading
elements of Markov transition probabilities from
parameter vector @
local pm,ixth;
ixth = rows(xth);
pm = zeros(ns,ns);
if ipm == 1; @ for ns =2 this option has parameters as p11 and p22 @
if izz == 1;
pm[1,1] = xth[1,1]^2/(1 +xth[1,1]^2);
pm[2,2] = xth[2,1]^2/(1 + xth[2,1]^2);
else;
pm[1,1] = xth[1,1];
pm[2,2] = xth[2,1];
endif;
pm[2,1] = 1 - pm[1,1];
pm[1,2] = 1 - pm[2,2];
elseif ipm == 2; @ general case has parameters pij for i = 1,...,n and
j = 1,...,n-1 @
pm[1:ns-1,.] = reshape(xth[1:ixth,1],ns-1,ns);
if izz == 1;
pm[ns,.] = ones(1,ns);
pm = pm^2;
pm = pm./(sumc(pm)');
else;
pm[ns,.] = (1 - sumc(pm))';
endif;
elseif ipm == 3; @ This section can be rewritten by user to impose zeros
and ones where desired @
if izz == 1;
pm[1,1] = xth[1,1]^2/(1 + xth[1,1]^2);
pm[1,3] = xth[2,1]^2/(1 + xth[2,1]^2 + xth[4,1]^2);
pm[2,1] = 1/(1 + xth[1,1]^2);
pm[2,2] = xth[3,1]^2/(1 + xth[3,1]^2);
pm[2,3] = xth[4,1]^2/(1 + xth[2,1]^2 + xth[4,1]^2);
pm[3,2] = 1/(1 + xth[3,1]^2);
pm[3,3] = 1/(1 + xth[2,1]^2 + xth[4,1]^2);
elseif izz == 2;
pm[1,1] = xth[1,1];
pm[1,3] = xth[2,1];
pm[2,1] = 1 - xth[1,1];
pm[2,2] = xth[3,1];
pm[2,3] = xth[4,1];
pm[3,2] = 1 - xth[3,1];
pm[3,3] = 1 - xth[2,1] - xth[4,1];
endif;
endif;
retp(pm);
endp;
/* ======================================================================= */
@ This section echos values of parameters @
"Order of autoregression";;pphi;
"Order of ARCH process";;karch;
"Number of primitive states";;ns;
"Number of lagged states that affect y";;ps;
"First observation used for estimation is";;nk;
if larch == 0; "no leverage effect"; else; "with leverage effect"; endif;
if tarch == 0;"Distribution is Normal"; elseif tarch ==1;
"distribution is t"; endif;
proc echoo(th); @ proc to echo starting values @;
local spar,alpha0,phi,sig0,a0,a1,g1,pm,b1,l1,cm,nu,eps,rss;
alpha0 = th[1,1];
spar = 2;
if pphi > 0;
phi = th[2:2+pphi-1,1];
spar = spar+1;
else;
phi = 0;
endif;
"Constant term in regression";; alpha0;
"Autoregressive coefficients in regression";; phi';
if barch == 0;
sig0 = abs(th[spar,1]);
"Initial variance";;sig0;
spar = spar + 1;
else;
"Initial variance not neeeded ";
endif;
a0 = abs(th[spar,1]);
a1 = abs(th[spar+1:spar+karch,1]);
"Constant term in ARCH process";;a0;
"Coefficients on lagged epsilon squared in ARCH process";;a1';
spar = spar+1+karch;
pm = matpm(th[spar:spar+(ns*(ns-1))-1-irs,1]);
"(Transposed) matrix of transition probabilities";;pm;"";
"The state with no adjustment to ARCH process is state 1, with transition";
"probability ";;pm[1,1];
spar = spar+(ns*(ns-1))-irs;
g1 = abs(th[spar:spar+ns-2,1]);
"Vector of variance factors for states 2 through";;ns;;g1';
spar = spar+ns-1;
if larch == 1;
l1 = abs(th[spar,1]);
"Coefficient on negative lagged change for asymmetric effect";;l1;
spar = spar+1;
endif;
if tarch == 1;
nu = 2 + abs(th[spar,1]);
"degree of freedom for t distribution is";; nu;
spar = spar + 1;
endif;
retp(a0);
endp;
/* ================================================================ */
@ This section calls main programs @
hp = pattern1;
#include procs;
x = startval;
@ The following lines are for convenience of analysis and should be removed
for final calculations
izz = 2;
kc = 2;
ks = 2; @
call echoo(x);
"";"Initial values:";; x';
"Initial value for negative log likelihood:";; ofn(x);
"";"Do you wish to continue (y or n)?";;
zzs = cons;
if zzs .$== "n";
end;
endif;
@ goto jump; @
/* ==================================================================== */
@ Set parameters to use Gauss numerical optimizer @
library optmum;
#include optmum.ext;
__btol = 1.e-08; @ This controls convergence criterion for coefficients@
__gtol = 1.e-08; @ This controls convergence criterion for gradient @
__algr = 1; @ This chooses BFGS optimization @
__miter = 80; @ This controls the maximum number of iterations @
__output = 1; @ This causes extra output to be displayed @
__covp = 0; @ This speeds up return from OPTMUM @
@ Next call the GAUSS numerical optimizer @
{x,f,g,h} =optmum(&ofn,startval);
"";"";"======================================================";
" FINAL ESTIMATES";
"";"Value of log likelihood:";;-f;
"";"Coefficients:";x';"";
call echoo(x);
"";"Gradient vector:";g';
vg=hessp(&ofn,x);
{va,ve} = eigrs2(vg);
va = sortc(va~ve',1);
if va[1,1] > 0;
"Standard errors:";
h=invpd(vg);
hh=sqrt(diag(h));
hh';
else;
"Hessian not positive definite; eigenvalues are";
va[.,1]';
"eigenvector associated with smallest eigenvalue is";
va[1,.];
endif;
/* ======================================================================= */
@ Print out complete analysis @
jump:
kc = 2;
ks = 2;
call ofn(x);
nxx = captst; @ Use nxx = captst for full output @
"Probabilities for primitive states";
"filtered probabilities";format /rd 1,0;
"Obs ";;"return ";;
t = 0;
do until t > ps;
i = 1;
do until i == ns;
"P(st-";;t;;"=";;i;;") ";;
i = i+1;
endo;
t = t+1;
endo;"";
format /rd 6,4;
skif = (skif*hp')*(eye(ps+1).*.id[.,1:ns-1]);
skif = seqa(nk,1,captst)~y[nk:capt,1]~skif;
skif[1:nxx,.];
"";"smoothed probabilities";
format /rd 1,0;
"Obs ";;"return ";;
i = 1;
do until i > ns;
"P(st = ";;i;;") ";;
i = i+1;
endo;
format /rd 6,4;
skis = skis*hp';
skis = seqa(nk,1,captst)~y[nk:capt,1]~skis[.,1:ns];
skis[1:nxx,.];
"";" Obs Residual Variance Prob of observing larger value";
outprob = 1 - outprob;
iitlow = 0;
iithigh = 0;
it = 1;
do until it > nxx;
nk+it-1;;" ";;resid[nk+it-1,1];;" ";;varfor[nk+it-1,1];;" ";;
if outprob[it+nk-1,1] > 0.975;
outprob[it+nk-1,1];;" *";
iitlow = iitlow+1;
elseif outprob[it+nk-1,1] < 0.025;
outprob[it+nk-1,1];;" *";
iithigh = iithigh+1;
else;
outprob[it+nk-1,1];
endif;
it = it+1;
endo;
"";"";"Number of observations below .025 level:";iitlow;
"Number of observations above .975 level:";iithigh;
/* hh = sortc(skis[1:captst,1]~outprob[nk:capt,1],2);
"";"Observations with lowest p values";
"";" Obs P-value";
format /m1 /ros 16,8;
hh[1:25,.];
"";"Observations with highest p values"; */
format /m1 /ros 16,8;
hh = skis[1:captst,1]~outprob[nk:capt,1];
hh = sortc(hh,2);
hh;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -