subroutine  backtests

pageselect quarterly

scalar h_min = 1
scalar h_max = 12

scalar model_max = 36

scalar n_grid = 1000  ' Grid for kernel pdf estimation

scalar T_VaR = 109  ' T_VaR=96 --> h=1 (1985q1:08q4) 

scalar pvalue_uc
scalar LR_uc
matrix(1,1) m1
matrix(1,1) m2
matrix(1,1) m3

matrix(9,model_max) zh1_percent_hits
matrix(9,model_max) zh1_p_value_Kupiec
matrix(9,model_max) zh1_p_value_Christoff
matrix(9,model_max) zh2_percent_hits
matrix(9,model_max) zh2_p_value_Kupiec
matrix(9,model_max) zh2_p_value_Christoff
matrix(9,model_max) zh3_percent_hits
matrix(9,model_max) zh3_p_value_Kupiec
matrix(9,model_max) zh3_p_value_Christoff
matrix(9,model_max) zh4_percent_hits
matrix(9,model_max) zh4_p_value_Kupiec
matrix(9,model_max) zh4_p_value_Christoff
matrix(9,model_max) zh8_percent_hits
matrix(9,model_max) zh8_p_value_Kupiec
matrix(9,model_max) zh8_p_value_Christoff
matrix(9,model_max) zh12_percent_hits
matrix(9,model_max) zh12_p_value_Kupiec
matrix(9,model_max) zh12_p_value_Christoff

matrix(9,model_max) zh1_p_value_vqr
matrix(9,model_max) zh2_p_value_vqr
matrix(9,model_max) zh3_p_value_vqr
matrix(9,model_max) zh4_p_value_vqr
matrix(9,model_max) zh8_p_value_vqr
matrix(9,model_max) zh12_p_value_vqr

matrix(9,model_max) zh1_stat_indep
matrix(9,model_max) zh2_stat_indep
matrix(9,model_max) zh3_stat_indep
matrix(9,model_max) zh4_stat_indep
matrix(9,model_max) zh8_stat_indep
matrix(9,model_max) zh12_stat_indep

'-------------------------------------------------------------------------------------------------
' *** Generation of inflation series, for each forecast horizon
for !h=h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then

scalar i_max = 165  '  up to 2012q1   
for !i=1 to i_max

pageselect Vintages_GDP_price_index

if (!i<>168) then
smpl @all
scalar calc_vint = !i + 2     ' Following Clark (2011, section 2.Data), we use the second available estimates of GDP/GNP and the GDP/GNP deflator as actuals in evaluating forecast accuracy. In the case of h-step ahead forecasts made for period t + h with vintage t data ending in period t - 1, the second available estimate is normally taken from the vintage t + h + 2 data set.

if (calc_vint=3) then
  scalar year = 1970
  scalar quarter_ = 4
endif

if (calc_vint=4) then
  scalar year = 1971
  scalar quarter_ = 1
endif

if (calc_vint>4) then
  scalar calc_years = @floor(calc_vint/4 )
  scalar calc_quarters = calc_vint - calc_years*4
  scalar year = 1970 + calc_years                            
  scalar quarter_ = 1 + calc_quarters
  if (quarter_>4) then
    scalar calc_aux = @floor(quarter_/4)
    scalar quarter_ = quarter_ - calc_aux*4
    scalar year = year + calc_aux
  endif
endif

for !s=year to year
for !q= quarter_ to quarter_
   series gdp_price_vintage = gdp_price_p!sq!q 
next
next

endif

if (!i=168) then
   smpl @all                                
   series gdp_price_vintage = gdp_price_p2012q2
endif

'--------------------------------------------------------------------------------------------

pageselect quarterly
smpl @all
series pi_a = 400*log(gdp_price_vintage/gdp_price_vintage(-1))
series pi_b = 400*log(gdp_price_vintage/gdp_price_vintage(-1))    '    100*log(gdp_price_vintage/gdp_price_vintage(-4))

smpl 1970q1+!i        1970q1+!i    '+!h-1    '+!h-1  

'if (!h>4) then
'    smpl 1947q1       1970q1+!i    '+!h-1    
'endif
    

if (!h<4) then
   series pi_aux_h!h = pi_a
   series pi_aux_vqr_h!h = log(1+pi_a/100)  -  log(1+pi_a(-!h)/100)
endif
if (!h=4) then
   series pi_aux_h!h = pi_b
   series pi_aux_vqr_h!h = log(1+pi_b/100) - log(1+pi_b(-!h)/100)
endif
if (!h>4) then
   series pi_aux_h!h = pi_b
   series pi_aux_vqr_h!h = log(1+pi_b/100) - log(1+pi_b(-!h)/100)
endif

next   ' !i


endif  ' h<5, =8, =12
next   ' !h


' end of loopings

'--------------------------------------------------------------------------------------------------------------------------------------------------------
' Forecast errors and MSE (models 1 to 16) to compute combination weights MSEW (Kascha & Ravazzolo, 2010)

smpl @all
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to 16   'models
   series forec_error_m!i_h!h =   (pi_aux_h!h - m!i_h!h_q50)
   series forec_error2_m!i_h!h = (pi_aux_h!h - m!i_h!h_q50) ^2
   series MSEWg1_m!i_h!h  ' Set of models 1-8
   series MSEWg2_m!i_h!h  ' Set of models 9-16
   series MSEWg3_m!i_h!h  ' Set of models 1-16
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

'--------------------------------------------------------------------------------------------------------------------------------------------------------
' Evaluate MSE along the Training Sample (TS)

matrix(model_max,h_max) fe_mean
matrix(model_max,h_max) fe_inv_mse
matrix(model_max,h_max) fe_MSEW
matrix(1,1) aux_msew

for !t = 1 to 109        ' 1985q1 to 2012q1
  smpl 1985q1-TS  1985q1+!t-1

for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then

for !i = 1 to 16
   scalar aux = @mean(forec_error_m!i_h!h)
   fe_mean(!i,!h) = aux
   scalar aux = @mean(forec_error2_m!i_h!h)
   fe_inv_mse(!i,!h) = 1 / aux   '  inverse of MSE and NOT root MSE (RMSE).
next  ' !i

endif  ' h<5, =8, =12
next  ' !h

' ---------Set of models 1-8 -----------------
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
    scalar sum_inv = 0
for !i = 1 to 8
    sum_inv = sum_inv + fe_inv_mse(!i,!h)
next  ' !i
for !i = 1 to 8
     fe_MSEW(!i,!h) = fe_inv_mse(!i,!h) / sum_inv 
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

smpl 1985q1+!t-1  1985q1+!t-1

for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to 8
        aux_msew(1,1) = fe_MSEW(!i,!h)
        mtos(aux_msew,MSEWg1_m!i_h!h)
next  ' i
endif  ' h<5, =8, =12
next  ' h

' ------------- Set of models 9-16 --------------
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
    scalar sum_inv = 0
for !i = 9 to 16
    sum_inv = sum_inv + fe_inv_mse(!i,!h)
next  ' !i
for !i = 9 to 16
     fe_MSEW(!i,!h) = fe_inv_mse(!i,!h) / sum_inv 
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

smpl 1985q1+!t-1  1985q1+!t-1

for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 9 to 16
        aux_msew(1,1) = fe_MSEW(!i,!h)
        mtos(aux_msew,MSEWg2_m!i_h!h)
next  ' i
endif  ' h<5, =8, =12
next  ' h

' -------------- Set of models 1-16 ---------------
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
    scalar sum_inv = 0
for !i = 1 to 16
    sum_inv = sum_inv + fe_inv_mse(!i,!h)
next  ' !i
for !i = 1 to 16
     fe_MSEW(!i,!h) = fe_inv_mse(!i,!h) / sum_inv 
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

smpl 1985q1+!t-1  1985q1+!t-1

for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to 16
        aux_msew(1,1) = fe_MSEW(!i,!h)
        mtos(aux_msew,MSEWg3_m!i_h!h)
next  ' i
endif  ' h<5, =8, =12
next  ' h



next ' t


'--------------------------------------------


endsub
call backtests
