' 1. JOINT THRESHOLD & PARAMETERS ESTIMATION 
' [Caner, Hansen (2004)]

'Define the threshold variable
smpl 1999M01 2009M04
genr gamma=libor3m

'Define the forecast horizon
scalar k=12

'Define trimming parameter
scalar trim=.15

vector(99) tau
vector(99) ssr
vector(100) ssr2
ssr2(1)=10000000
vector(117) e_ch1_v
matrix(117,99) resids

'Generate 100 threshold values
smpl 1999M01 2009M04
for !j=1 to 99
tau(!j)=@quantile(gamma,trim)+!j*(@quantile(gamma,(1-trim))-@quantile(gamma,trim))/100

'Generate the corresponding dummy variables
smpl 1999M01 2009M04
genr d1=0
smpl if gamma(-1) > tau(!j)
genr d1=1

'Estimate the model conditionally on each gamma
smpl 1999M01 2009M03
equation cls.tsls libor3m c cpi*d1 gap*d1  fxgap*d1 cpi*(1-d1) gap*(1-d1) fxgap*(1-d1) @ cpi(-1 to -6, -9, -12) gap(-1 to -6, -9, -12) fxgap(-1 to -6, -9, -12)

'If necessary, check the behaviour of the residuals
'cls.makeresids e_ch1
'stom(e_ch1,e_ch1_v)
'colplace(resids,e_ch1_v,!j)

'Find tau_star
smpl 1999M01 2009M03
ssr(!j)=cls.@ssr
if ssr(!j) < ssr2(!j) then
genr tau_star=tau(!j)
ssr2(!j+1)=ssr(!j)
else
ssr2(!j+1)=ssr2(!j)

endif
next

' Re-estimate the model conditional on tau_star separately for each regime
smpl if gamma(-1) > tau_star
equation ccls1.gmm libor3m c cpi(+k) gap fxgap @ cpi(-1 to -6, -9, -12) gap(-1 to -6, -9, -12) fxgap(-1 to -6, -9, -12)
freeze(mode=overwrite, tableccls1) ccls1.gmm
scalar obs1=@obs(gamma(-1))

smpl if gamma(-1) < tau_star
equation ccls2.gmm libor3m c cpi(+k) gap fxgap @ cpi(-1 to -6, -9, -12) gap(-1 to -6, -9, -12) fxgap(-1 to -6, -9, -12)
freeze(mode=overwrite, tableccls2) ccls2.gmm
scalar obs2=@obs(gamma(-1))

'Prepare a matrix to store the results of interest
matrix(8,8) ci

'Store the point estimates of the slope parameters
ci(1,1)=@val(tableccls1(14,2))
ci(2,1)=@val(tableccls1(15,2))
ci(3,1)=@val(tableccls1(16,2))
ci(4,1)=@val(tableccls1(17,2))

ci(1,5)=@val(tableccls2(14,2))
ci(2,5)=@val(tableccls2(15,2))
ci(3,5)=@val(tableccls2(16,2))
ci(4,5)=@val(tableccls2(17,2))

'Store the std error of the slope parameters
ci(1,4)=@val(tableccls1(14,3))
ci(2,4)=@val(tableccls1(15,3))
ci(3,4)=@val(tableccls1(16,3))
ci(4,4)=@val(tableccls1(17,3))

ci(1,8)=@val(tableccls2(14,3))
ci(2,8)=@val(tableccls2(15,3))
ci(3,8)=@val(tableccls2(16,3))
ci(4,8)=@val(tableccls2(17,3))

'Store the estimate threshold*
ci(5,1)=tau_star(1)

'Compute and store the implicit inflation target pi*
ci(7,1)=(0.6-@val(tableccls1(14,2)))/(@val(tableccls1(15,2))-1)
ci(7,5)=(0.6-@val(tableccls2(14,2)))/(@val(tableccls2(15,2))-1)

'Compute and store the p-values of the J-Statistic
scalar overid1=ccls1.@regobs*ccls1.@jstat
ci(8,1)=1-@cchisq(overid1,23)
scalar overid2=ccls2.@regobs*ccls2.@jstat
ci(8,5)=1-@cchisq(overid2,23)

'Compute the number of observations in each regime
ci(6,1)=obs1
ci(6,5)=obs2


' 2. CONFIDENCE INTERVALS
' [Hansen (2000), Caner & Hansen (2004)]

' Compute visually confidence interval for tau_star
smpl 1999M01 2007M03
series ssrs
series taus
smpl 1999M01 2007M03
mtos(ssr,ssrs)
mtos(tau,taus)
series ssrmi=@min(ssrs)
series lrss=124*(ssrs-ssrmi)/ssrmi
smpl 1999M01 2007M03
xyline taus lrss

'Compute numerically 95% confidence interval for tau_star
'critical value=7.35 {Hansen, 1997]
'Beware: only valid is the interval is not disjoint
smpl 1999M01 2007M03 if lrss<7.35
ci(5,2)=@min(taus)
ci(5,3)=@max(taus)

'Compute 80% confidence interval for the slope parameters
'critical value=4.5 {Hansen, 1997]

vector(obsci) cilc1
vector(obsci) cilcpi1
vector(obsci) cilgap1
vector(obsci) cilfxgap1
vector(obsci) ciuc1
vector(obsci) ciucpi1
vector(obsci) ciugap1
vector(obsci) ciufxgap1
vector(obsci) cilc2
vector(obsci) cilcpi2
vector(obsci) cilgap2
vector(obsci) cilfxgap2
vector(obsci) ciuc2
vector(obsci) ciucpi2
vector(obsci) ciugap2
vector(obsci) ciufxgap2

'Select the gamma that lie within the 80% confidence interval
smpl 1999M01 2007M03 if lrss<4.5
scalar obsci=@obs(taus)
vector(obsci) tauci

'Re-estimate the models for the gamma that lie within the 80% confidence interval -> Regime 1
stom(taus,tauci)
for !k=1 to obsci
smpl 1999M01 2009M03
smpl if gamma(-1) > tauci(!k)
equation acils1.gmm libor3m c cpi(+k) gap fxgap @ cpi(-1 to -6, -9, -12) gap(-1 to -6, -9, -12) fxgap(-1 to -6, -9, -12)
freeze(mode=overwrite, tableci1) acils1.gmm

'Compute the 95% confidence interval for each parameter in each model
vector cilc1(!k)=@val(tableci1(14,2))-1.96*@val(tableci1(14,3))
vector cilcpi1(!k)=@val(tableci1(15,2))-1.96*@val(tableci1(15,3))
vector cilgap1(!k)=@val(tableci1(16,2))-1.96*@val(tableci1(16,3))
vector cilfxgap1(!k)=@val(tableci1(17,2))-1.96*@val(tableci1(17,3))

vector ciuc1(!k)=@val(tableci1(14,2))+1.96*@val(tableci1(14,3))
vector ciucpi1(!k)=@val(tableci1(15,2))+1.96*@val(tableci1(15,3))
vector ciugap1(!k)=@val(tableci1(16,2))+1.96*@val(tableci1(16,3))
vector ciufxgap1(!k)=@val(tableci1(17,2))+1.96*@val(tableci1(17,3))

'Compute the joint confidence interval
ci(1,2)=@min(cilc1)
ci(2,2)=@min(cilcpi1)
ci(3,2)=@min(cilgap1)
ci(4,2)=@min(cilfxgap1)

ci(1,3)=@max(ciuc1)
ci(2,3)=@max(ciucpi1)
ci(3,3)=@max(ciugap1)
ci(4,3)=@max(ciufxgap1)

'Do the same for the second regime
smpl 1999M01 2009M03
smpl if gamma(-1) < tauci(!k)
equation acils2.gmm libor3m c cpi(+k) gap fxgap @ cpi(-1 to -6, -9, -12) gap(-1 to -6, -9, -12) fxgap(-1 to -6, -9, -12)
freeze(mode=overwrite, tableci2) acils2.gmm

vector cilc2(!k)=@val(tableci2(14,2))-1.96*@val(tableci2(14,3))
vector cilcpi2(!k)=@val(tableci2(15,2))-1.96*@val(tableci2(15,3))
vector cilgap2(!k)=@val(tableci2(16,2))-1.96*@val(tableci2(16,3))
vector cilfxgap2(!k)=@val(tableci2(17,2))-1.96*@val(tableci2(17,3))
ci(1,6)=@min(cilc2)
ci(2,6)=@min(cilcpi2)
ci(3,6)=@min(cilgap2)
ci(4,6)=@min(cilfxgap2)

vector ciuc2(!k)=@val(tableci2(14,2))+1.96*@val(tableci2(14,3))
vector ciucpi2(!k)=@val(tableci2(15,2))+1.96*@val(tableci2(15,3))
vector ciugap2(!k)=@val(tableci2(16,2))+1.96*@val(tableci2(16,3))
vector ciufxgap2(!k)=@val(tableci2(17,2))+1.96*@val(tableci2(17,3))
ci(1,7)=@max(ciuc2)
ci(2,7)=@max(ciucpi2)
ci(3,7)=@max(ciugap2)
ci(4,7)=@max(ciufxgap2)

next