OPEN DATA "\\myfiles\2018 nonlinear paper revision\RATS_real time data for paper(shadow rate)_USUK.xls"
CALENDAR 1979 1 4
ALL 2015:04
DATA(FORMAT=XLS,ORG=COLUMNS) 1979:01 2015:04 uk_ex  uk_int	us_int $
	uk_qdy us_qdy uk_inf us_inf uk_realex us_qds	uk_qds	uk_qdh	us_qdh

print 2005:1 2006:1 us_inf


****************************************************
*data transform

set y = uk_ex
set diffinf =us_inf - uk_inf
set diffqdy = us_qdy-uk_qdy
set diffint = us_int-uk_int


set diffws = us_qds-uk_qds
set diffwh = us_qdh-uk_qdh
set realex = uk_realex

****************************************************

*Linearity test on STR models

set y = uk_ex
set a1 = diffinf{1}
set a2 = diffint{2}
set a4 = diffqdy{1}
set a5 = diffws{1}
set a6 = diffwh{1}
set a7 = realex{1}


linreg y 1980:1 2015:4
# constant diffinf{1} diffint{2} diffqdy{1} diffwh{1} diffws{1} realex{1}
@regSTRTest(threshold=a1, d=0) y
@regSTRTest(threshold=a2, d=0) y
@regSTRTest(threshold=a4, d=0) y
@regSTRTest(threshold=a5, d=0) y
@regSTRTest(threshold=a6, d=0) y
@regSTRTest(threshold=a7, d=0) y
@regSTRTest(threshold=volatility, d=0) y


****************************************************************************
***Estimation of the nonlinear model
******The ESTR model***
**
****************************************************************************
*
equation standard y
# constant diffinf{1} diffqdy{1} diffint{2} diffws{1} diffwh{1} realex{1}
equation transit y
# constant diffinf{1} diffqdy{1} diffint{2} diffws{1} diffws{1} realex{1}
*
dec vector phi1 phi2
frml(equation=standard,vector=phi1) phi1f
frml(equation=transit,vector=phi2) phi2f
nonlin(parmset=regparms) phi1 phi2
*
nonlin(parmset=starparms) gamma c
*
frml festar   = 1.0-exp(-gamma*(diffint{2}-c)^2)
frml estar y = g=festar,phi1f+g*phi2f


stats diffint
compute gamma=1.0/%variance
compute c = %mean

********************************************************************
* Estimate the model with the ESTAR parameters fixed.
*
nlls(parmset=regparms,frml=estar) y
*
* Then with the STAR parameters freed
*
nlls(parmset=regparms+starparms,frml=estar) y / resids
********************************************************************
*Grid search on c and gamma

compute agrid=%seqa(1, 1, 10)
dec real ggrid

stats(fractiles) diffint
compute grid=%seqa(%fract25,(%fract75-%fract25)/49.0,50)
dec real cgrid

compute minrss=1.e+100
nonlin(parmset=starparms) gamma=ggrid c=cgrid

dofor ggrid =agrid
    dofor cgrid = grid
      nlls(frml=estar,parmset=regparms+starparms,noprint) y
       if %rss<minrss
           compute minrss=%rss,cbest=cgrid, gbest=ggrid
   disp c gamma %rss
    end dofor cgrid
end dofor ggrid
dis 'c = ' cbest
dis 'gamma = ' gbest


********************************************************************
* Estimates a final model starting at the best value for gamma and c

compute cstar = cbest
nonlin(parmset=starparms) gamma cstar
nlls(frml=estar,parmset=regparms+starparms,print) y


****************************************************************************
***Estimation of the nonlinear model
******The LSTR model***
**
****************************************************************************
equation standard y
# constant diffinf{1} diffqdy{1} diffint{2} diffws{1} diffwh{1} realex{1}
equation transit y
# constant diffinf{1} diffqdy{1} diffint{2} diffws{1} diffwh{1} realex{1}
*
dec vector phi1 phi2
frml(equation=standard,vector=phi1) phi1f
frml(equation=transit,vector=phi2) phi2f
nonlin(parmset=regparms) phi1 phi2
*
nonlin(parmset=starparms) gamma c
*
frml glstar = %logistic(gamma*(diffint{2}-c),1.0)
frml star y = g=glstar,phi1f+g*phi2f


stats diffint
compute gamma=2.0/sqrt(%variance)
compute c = %mean

********************************************************************
* Estimate the model with the LSTAR parameters fixed.
*
nlls(parmset=regparms,frml=star) y
*
* Then with the STAR parameters freed
*
nlls(parmset=regparms+starparms,frml=star) y / resids
   set u = %resids
   print 1976:3 1976:4 u

test
# 7 8 9 10 11 12
# 0 0 0 0 0 0


********************************************************************
*Grid search on c and gamma

compute agrid=%seqa(1, 1, 10)
dec real ggrid

stats(fractiles)  diffint
compute grid=%seqa(%fract25,(%fract75-%fract25)/49.0, 50)
dec real cgrid

compute minrss=1.e+100
nonlin(parmset=starparms) gamma=ggrid c=cgrid

dofor ggrid =agrid
    dofor cgrid = grid
      nlls(frml=star,parmset=regparms+starparms,noprint) y
       if %rss<minrss
           compute minrss=%rss,cbest=cgrid, gbest=ggrid
   disp c gamma %rss
    end dofor cgrid
end dofor ggrid
dis 'c = ' cbest
dis 'gamma = ' gbest

********************************************************************
* Estimates a final model starting at the best value for gamma and c


compute cstar = cbest
compute gamma = gbest
nonlin(parmset=starparms) gamma cstar
nlls(frml=star,parmset=regparms+starparms,print) y / resids
   set u = %resids
   print 1976:3 1976:4 u




********************************************************************
***Non-Linear forecasting
***** MSPE
****************************************************************************

group  starmodel star>>sforecast
set meanf 1989:1 2008:4 = 0.0
compute nreps=10000
do reps=1,nreps
   simulate(model=starmodel,from=1989:1,to=2008:4)
   set meanf 1989:1 2008:4 = meanf+sforecast
end do reps
set meanf 1989:1 2008:4 = meanf/nreps


forecast(model=starmodel, from=1989:1,to=2008:4, print)


@UforeErrors y sforecast 1989:1 2008:4
@UforeErrors  y meanf 1989:1 2008:4


***********************************************************************
***Linear Forecasting (benchmark: RW)
***** MSPE, CW
****************************************************************************
*UK

***** fixed window rolling regression( Model 1)******

clear uk_fex coef1 coef2 coef3 coef4 coef5 coef6 coef7 coef8 coef9 coef10 coef11
do regend = 1988:4, 2008:4
    linreg(noprint) uk_ex regend-39 regend
    # Constant diffinf{1} diffhpy{1} diffws{1} uk_realex{1} diffint{2}
    compute coef1(regend) = %beta(1)
    compute coef2(regend) = %beta(2)
    compute coef3(regend) = %beta(3)
    compute coef4(regend) = %beta(4)
    compute coef5(regend) = %beta(5)
    compute coef6(regend) = %beta(6)
    uforecast uk_fex1 regend+1 regend+1
end do

@UforeErrors(title="usuk(1) forecasts")  uk_ex uk_fex1 1989:1 2008:4


*******CW Statistics (RW)********

SET FORE1 1989:01 2008:04 = 2*uk_ex*uk_fex1
linreg(robusterrors,lags=4,lwindow=newey) fore1 1989:1 2008:4
# constant


********************************************************************
***Non-Linear forecasting (benchmark: linear Model & RW)
*****CW
****************************************************************************

SET FORE1 1989:01 2008:04 = 2*(uk_fex1-y)*(sforecast-uk_fex1)
linreg fore1 1989:1 2008:4
# constant

SET FORE1 1989:01 2008:04 = 2*(uk_fex1-y)*(meanf-uk_fex1)
linreg fore1 1989:1 2008:4
# constant

SET FORE1 1989:01 2008:04 = 2*y*sforecast
linreg fore1 1989:1 2008:4
# constant