Here is the code which is an adjusted version of your code:
Code: Select all
*
* All of this was out of Gray's paper, "Modeling the conditional
* distribution of interest rates as a regime-switching process", J. of
* Financial Economics 42, 1996, pp 27-62.
*
* Revision Schedule:
* 01/2005 Written by Tom Doan, Estima
* 09/2005 Comments added. Changed to use new %MSSMOOTH procedure
* 08/2011 Changed to use the new MSSETUP procedures.
*
OPEN DATA "D:sf.xlsx"
CALENDAR(D) 2002:1:3
DATA(FORMAT=XLSX,ORG=COLUMNS) 2002:01:03 *** data rate
set rate = rate*100.0
spgraph(vfields=2,window="Figure 3")
graph(header=" Exchnage Rate")
# data
graph(header="Daily Exchange Rate Return")
# rate
spgraph(done)
*
* This makes extensive use of the Markov switching functions and
* procedures from MSSetup.src
*
@MSSetup(states=2)
*
* The P and Q values in Gray's notation are handled in MSSetup as an M-1
* x M matrix (called P here), with P(i,j) giving the probabilities of
* moving from state j to state i. The Mth row is implicit, as it's just
* the 1-sum of rest of the column.
*
*******************************************************************************
*
* Table 1 information
*
stats rate
cmom(corr,print)
# rate rate{1}
*******************************************************************************
*
* Table 2 estimates
*
* Single regime is linear regression
*
linreg(robust) rate
# constant rate{1}
*
* We're going to use these later for initial conditions
*
compute olsvar =%seesq
compute olsbeta=%beta
*
* Constant variance regime switching model
* This uses Gray's notation for all the variables other than "P"
*
nonlin a01 a02 a11 a12 b01 b02 p
*
* Initial guess values for MS models are always a bit tricky because the
* model isn't globally identified. This is starting with OLS results for
* state 1 and 0 mean coefficients for state 2.
*
compute a01=olsbeta(1),a11=olsbeta(2),b01=sqrt(olsvar)
compute a02=a12=0.0,b02=b01
*
* The p matrix is initialized to make both states fairly persistent.
* (The 2,2 transition is 1-.1=.9).
*
compute p=||.8,.1||
*
* SimpleRegimeF returns the 2-vector of densities in the two regimes. It
* is specific to the data series used (dependent variable <<rate>>) and
* the precise parameterization, and so would need to be altered slightly
* for a different application. Note that this needs to return the actual
* density, not its log.
*
function SimpleRegimeF t
type vector SimpleRegimeF
type integer t
*
compute SimpleRegimeF=||%density((rate(t)-a01-a11*rate(t-1))/b01)/b01,$
%density((rate(t)-a02-a12*rate(t-1))/b02)/b02||
end
*
* This is the standard log likelihood FRML for a Markov switching model.
*
frml logl = f=SimpleRegimeF(t),fpt=%MSProb(t,f),log(fpt)
@MSFilterInit
maximize(start=(pstar=%MSInit()),method=bfgs,iters=100,pmethod=simplex,piters=5) logl 2 ***
*
* It's not clear what residuals are used for the diagonostics for the
* regime switching model. The standardized residuals weighted by the ex
* ante probabilities of the two states seem to give LB's which are
* fairly close to those in the paper.
*
set ustd = %dot(pt_t1,||(rate-a01-a11*rate{1})/b01,(rate-a02-a12*rate{1})/b02||)
graph
# ustd
set usqr = ustd^2
@regcorrs(report,number=15) usqr
*
* Compute the smoothed probabilities and the ex ante probabilities
*
@MSSmoothed %regstart() %regend() psmooth
set exante %regstart() %regend() = pt_t1(t)(1)
set smooth %regstart() %regend() = psmooth(t)(1)
*
* And the conditional standard deviation
*
set condstddev = sqrt(exante*b01^2+(1-exante)*b02^2)
*
spgraph(vfields=2,window="Figure 4")
graph(max=1.0,header="Regime Probabilities") 2
# exante
# smooth
graph(header="Conditional Standard Deviation")
# condstddev
spgraph(done)
*******************************************************************************
*
* Table 3 estimates
*
* Single regime GARCH model
*
garch(p=1,q=1,reg,resids=u,hseries=h) / rate
# constant rate{1}
compute onestate=%beta
*
set usqr = u^2/h
@regcorrs(number=15,report) usqr
set usqr1 = u/h
@regcorrs(number=15,report) usqr1
*
* Regime-switching GARCH
*
nonlin a01 a02 a11 a12 b01 b11 b21 b02 b12 b22 p
compute a01=%beta(1),a11=%beta(2),b01=olsvar,b11=0.0,b21=0.0
compute a02=%beta(1),a12=%beta(2),b02=%beta(3),b12=%beta(4),b22=%beta(5)
compute p=||.95,.05||
*
set uu = olsvar
set h = olsvar
*
* RegimeGARCHF returns the 2 vector of densities in the two regimes.
* Again, many lines in this are specific to the problem at hand.
*
function RegimeGARCHF t
type vector RegimeGARCHF
type integer t
local real hh1 hh2 mu1 mu2 mu
*
* Compute state dependent variances
*
compute hh1=b01+b11*uu(t-1)+b21*h(t-1)
compute hh2=b02+b12*uu(t-1)+b22*h(t-1)
*
* Compute state dependent means
*
compute mu1=a01+a11*rate(t-1)
compute mu2=a02+a12*rate(t-1)
*
* Compute the return vector (densities in the two states)
*
compute RegimeGARCHF=||%density((rate(t)-mu1)/sqrt(hh1))/sqrt(hh1),$
%density((rate(t)-mu2)/sqrt(hh2))/sqrt(hh2)||
*
* Compute the values of uu (squared residual) and h (variance) to be
* used for the period following.
*
compute mu=mu1*pstar(1)+mu2*pstar(2)
compute uu(t)=(rate(t)-mu)^2
compute h(t)=pstar(1)*(mu1^2+hh1)+pstar(2)*(mu2^2+hh2)-mu^2
end
*
frml logl = f=RegimeGARCHF(t),fpt=%MSProb(t,f),log(fpt)
*
* This combination is able to locate Gray's results - the first MAXIMIZE
* works with the "p" matrix fixed at a fairly high level of persistence,
* and tries to get estimates which separate the two states. The second
* then adds the p matrix to the parmset.
*
nonlin a01 a02 a11 a12 b01 b11 b21 b02 b12 b22
maximize(start=(pstar=%MSInit()),method=bfgs,iters=100,pmethod=simplex,piters=50) logl 2 ***
nonlin a01 a02 a11 a12 b01 b11 b21 b02 b12 b22 p
maximize(start=(pstar=%MSInit()),method=bfgs,iters=100) logl 2 ***
set usqr2 2 * = uu(t)/h(t)
@regcorrs(title="Squared Standardized Residuals from MS GARCH",$
number=15,report) usqr2
graph
# usqr2
@regcorrs(report,number=15) usqr2
@westchotest(number=5) usqr2
@westchotest(number=10) usqr2
@westchotest(number=15) usqr2
set usqr3 = usqr2^.5
@regcorrs(title="Standardized Residuals from MS GARCH",$
number=15,report) usqr3
graph
# usqr3
*
* Compute the smoothed probabilities and the ex ante probabilities
*
@MSSmoothed %regstart() %regend() psmooth
set exante %regstart() %regend() = pt_t1(t)(1)
set smooth %regstart() %regend() = psmooth(t)(1)
*
* And the conditional standard deviation
*
set condstddev = sqrt(h)
*
spgraph(vfields=2,window="Figure 5")
graph(max=1.0,header="Regime Probabilities") 2
# exante
# smooth
graph(header="Conditional Standard Deviation")
# condstddev
spgraph(done)
*
*
* However, this set of initial guess values (which are actually the
* result of a typo in trying to reproduce the article's numbers) find a
* local max with quite a bit higher likelihood, and dramatically
* different behavior. How one might find this systematically is unclear.
*
compute a01=.1407,a02=-.0011,a11=-.0141,a12=.0006,b01=.0004,b11=.4609,b21=.1977,b02=.0099,b12=.1655,b22=.2685,p=||.9739,1-.9896||
maximize(start=(pstar=%MSInit()),method=bfgs,iters=100,pmethod=simplex,piters=50) logl 2 ***
*******************************************************************************