' Regime switching for Beta switching (Example 1.9,  Section 1.3). Regimes 0 = low volatility, 1 = high volatility

series y = fin       ' <<<<<<<<<<< change for a proper series name
series x = world   ' <<<<<<<<<<< Average mkt index. Change for a proper series

series  d1 = 0     ' dummy variable to control for the first observation

smpl @first + 2  @first + 2
d1 = 1

smpl @first + 2  @last                 ' Sample period. Two first observations lost in data adjustments

' Initial values
equation init.ls y c x     ' Create temprary init equation to compute initial estimates

' Define the log-likelihood object

logl msw_beta 

' Initial values for the parameter estimates

coef(1) alpha0 = init.c(1)            ' mean of the process in "low" volatility regime
coef(1) alpha1 = init.c(1)            ' mean of the process in "high" volatitlity regime
coef(1) beta0(1) = init.c(2)          ' Beta parameter estimate, "low" volatility regime
coef(1) beta1(1) = init.c(2)          ' Beta parameter estimate, l"high" volatility regime
coef(1) sigma0p2 = (0.8*init.@se)^2  ' residual variance in the "low" volatility regime
coef(1) sigma1p2 = (1.2*init.@se)^2  ' residual variance in the  "high volatility regime

' Transition probabilities 
'coef(1) p11 = 0.8    ' P[S(t) = 1 | S(t-1) = 1]
'coef(1) p00 = 0.8    ' P[S(t) = 0 | S(t-1) = 0]

'Reparametrize transition probabilities in terms of a00 and a11 such that e.g. p00 = exp(a00) / (1 + exp(a00))
' This quarantees that the probabilities in the numerical agorithm will stay within 0 and 1.

coef(1) a00 = log(0.8/0.2)    ' p00 = P[S(t) = 0 | S(t-1) = 0] = exp(a00)/(1 + exp(a00)), P[S(t) = 1 | S(t-1) = 0] = 1 - P[S(t) = 0 | S(t-1) = 0]
coef(1) a11= log(0.8/0.2)     ' p11 = P[S(t) = 1 | S(t-1) = 1] = exp(a11)/(1 + exp(a11)), P[S(t) = 0 | S(t-1) = 1] = 1 - P[S(t) = 1 | S(t-1) = 1]


' Define the loglikelihood function name (loglrsw)

msw_beta.append @logl loglmsw_beta

' Specify here your mean models in the regimes
' Standardized mean models residuals to be used in the normal distributions. CHANGE these according to your specification

msw_beta.append sres00 = ( y - alpha0(1) - beta0(1)*x )/@sqr(sigma0p2(1))     ' From t-1 regime 0 to t regime 0
msw_beta.append sres01 = ( y - alpha1(1) - beta1(1)*x )/@sqr(sigma1p2(1))     ' From t-1 regime 0 to t regime 1
msw_beta.append sres10 = ( y - alpha0(1) - beta0(1)*x )/@sqr(sigma0p2(1))     ' From t-1 regime 1 to t regime 0
msw_beta.append sres11 = ( y - alpha1(1) - beta1(1)*x )/@sqr(sigma1p2(1))     ' From t-1 regime 1 to t regime 1

'***************** Do NOT chanee the specifications below this line except the logl object name (msw) if you chanaged it above. ************'

' Probabilities P[S(t-1) = 0 | F(t-1)] for calculating the log-likelihood.  Do not change these (except the logl object name if you chenged it above)!
' Remark: For the first observation P[S(t-1) = | F(t-1)] = (1 - p11) / (2 - p00 - p11)
' p1tm1 = P[S(t-1) = 0 | F(t-1)] = ( f(y(t-1) | S(t-1) = 0, S(t-2) = 0, F(t-2)) * P[S(t-1)=0|S(t-2)=0]*p0tm1(-1)
'                                   + f(y(t-1) | S(t-1) = 0, S(t-2) = 1, F(t-2))*P[S(t-1)=0|S(t-2)=1]*p1tm1(-1) )/ f(y(t-1)| F(t-2))

msw_beta.append p0tm1 =  @recode( d1=1, ( 1 - exp(a00(1)) /(1+exp(a00(1))) )/( 2 - exp(a11(1))/(1+exp(a11(1))) - exp(a00(1))/(1+exp(a00(1))) ), _
                   ( f00(-1)*exp(a00(1))/(1+exp(a00(1)))*p0tm1(-1) + f01(-1)*(1-exp(a11(1))/(1+exp(a11(1))))*p1tm1(-1) )/ fyt(-1) ) 
msw_beta.append p1tm1 =  1 - p0tm1

' Normal densities

msw_beta.append f00 = @dnorm(sres00)/@sqr(sigma0p2(1))  ' Density function when S(t-1) = 0 and S(t) = 0
msw_beta.append f01 = @dnorm(sres01)/@sqr(sigma0p2(1))  ' Density function when S(t-1) = 0 and S(t) = 1
msw_beta.append f10 = @dnorm(sres10)/@sqr(sigma1p2(1))  ' Density function when S(t-1) = 1 and S(t) = 0
msw_beta.append f11 = @dnorm(sres11)/@sqr(sigma1p2(1))  ' Density funciton when S(t-1) = 1 and S(t) = 1

' Likelihood function

msw_beta.append fyt =  f00*exp(a00(1))/(1 + exp(a00(1)))*p0tm1 + f01*( 1-exp(a11(1))/(1+exp(a11(1))) ) * p1tm1 _
                         + f10*( 1 - exp(a00(1))/(1 + exp(a00(1))) )*p0tm1 + f11*exp(a11(1))/(1 + exp(a11(1)))*p1tm1

' Log-likelihood

msw_beta.append loglmsw_beta = log(fyt)

' Estimate by maximum likelihood

msw_beta.ml(showopts, m=1000, c = 1e-5)
msw_beta.updatecoefs

' Smoothed probabilities Kim and Nelson (1999) pp. 68--69. (Kim's Smoothing Algorithm)

' P[S(t) =  j, S(t+1) = k] = P[S(t+1) = k | F(T)] * P[S(t) = j | F(t)] * P[S(t+1) = k | S(t) = j] / P[S(t+1) = k | F(t)]
' where P[S(t+1) = k | F(t) ] = Sum(i=0, 1) P[S(t+1) = k, S(t) = i | F(t)]  = Sum(i=0,1) P[S(t+1)  = k | S(t) = i] P[S(t) = i | F(t)]

series smoothpt1 = p1tm1

' Transition probabilities

scalar p00 = exp(a00(1))/(1 + exp(a00(1)))  ' P[S(t) = 0 | S(t-1) = 0]
scalar p11 = exp(a11(1))/(1 + exp(a11(1)))  ' P[S(t) = 1 | S(t-1) = 1]


smpl @all

' Number of observation in the sample

!n = @obssmpl

' Smoothed probabilities 

for !i = 1 to !n-1
    smpl @last-!i @last-!i 

    smoothpt1 =  (1-smoothpt1(1)) * p1tm1 * (1 - p11)/(p00*p0tm1 + (1-p11)*p1tm1) _
                   + smoothpt1(1) * p1tm1 * p11/((1-p00)*p0tm1 + p11*p1tm1)
next

smpl @all

series smoothpt0 = 1 - smoothpt1
  
' Create Figures and show output

line(d) fin world
line(d) p1tm1 smoothpt1
show msw_beta.output


'Cleanup 

delete y x d1
delete f00 f10 f01 f11 fyt
delete sres00 sres01 sres10 sres11
delete alpha0 alpha1 beta0 beta1 a00 a11 sigma0p2 sigma1p2