subroutine  backtests
' ------------------------------------------------------------
'  ****  Out-of-sample PDFs generation ****
delete y_aux
matrix(1,1) fan_chart_aux
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to model_max   

matrix(99,model_max) fan_chart_h!h
matrix(99,h_max)  fan_chart_m!i

smpl 2009q1+!h 2009q1+!h

for !tau = 1 to 99
   fan_chart_h!h(!tau,!i) = NA
   fan_chart_m!i(!tau,!h) = NA
   if (!i=1) then
       stom(m1_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=2) then
       stom(m2_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=3) then
       stom(m3_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=4) then
       stom(m4_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=5) then
       stom(m5_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=6) then
       stom(m6_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=7) then
       stom(m7_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=8) then
       stom(m8_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=9) then
       stom(m9_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=10) then
       stom(m10_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=11) then
       stom(m11_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=12) then
       stom(m12_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=13) then
       stom(m13_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=14) then
       stom(m14_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=15) then
       stom(m15_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=16) then
       stom(m16_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=17) then
       stom(m17_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=18) then
       stom(m18_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=19) then
       stom(m19_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=20) then
       stom(m20_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=21) then
       stom(m21_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=22) then
       stom(m22_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=23) then
       stom(m23_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=24) then
       stom(m24_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=25) then
       stom(m25_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=26) then
       stom(m26_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=27) then
       stom(m27_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=28) then
       stom(m28_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=29) then
       stom(m29_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=30) then
       stom(m30_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=31) then
       stom(m31_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=32) then
       stom(m32_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=33) then
       stom(m33_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=34) then
       stom(m34_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=35) then
       stom(m35_h!h_q!tau,fan_chart_aux)   
   endif
   if (!i=36) then
       stom(m36_h!h_q!tau,fan_chart_aux)   
   endif


   fan_chart_h!h(!tau,!i) = fan_chart_aux(1,1)
   fan_chart_m!i(!tau,!h) = fan_chart_aux(1,1)

next  ' tau

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

' end of loopings


'---------------------------------------------------------------------------
' Forecast errors (RMSFE relative to SPF forecasts)
smpl @all
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to model_max   
   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
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

matrix(model_max,h_max) fe_mean
matrix(model_max,h_max) fe_mean2
matrix(model_max,h_max) fe_rmse
matrix(model_max,h_max) fe_rmse_relative_spf

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

smpl 1985q1  2012q1      

for !i = 1 to model_max
   scalar aux = @mean(forec_error_m!i_h!h)
   fe_mean(!i,!h) = aux
   scalar aux = @mean(forec_error2_m!i_h!h)
   fe_mean2(!i,!h) = aux
   fe_rmse(!i,!h) = aux^0.5
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to model_max
   fe_rmse_relative_spf(!i,!h) = fe_rmse(!i,!h) /  fe_rmse(31,!h)   ' comparison with model 31
next  ' !i
endif  ' h<5, =8, =12
next  ' !h

' end of loopings


'----------------------------------------------------
' *** Berkowitz (JBES, 2001) test ***
'delete ut
'delete zt

matrix(model_max,h_max)   z_p_value_berkowitz

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

for !i = 1 to model_max

smpl @all
series  y =  pi_aux_h!h                   

series ut
scalar tau = 0
scalar nobs = T_VaR

for !tau = 1 to 99   
    series y_quant!tau = m!i_h!h_q!tau
next

matrix (1,1) y_aux
matrix (1,1) y_quant_aux
matrix(1,1) tab_aux

for !t = 1 to nobs
             smpl 1985q1+!t-1   1985q1+!t-1   

             tab_aux(1,1) = 0.5    
             mtos(tab_aux,ut)
             scalar achou = 0
                for !tau = 1 to 99    
                         stom(y,y_aux)
                         stom(y_quant!tau,y_quant_aux)
                         if (y_aux<y_quant_aux) and (achou=0)  then
                                 achou = 1                                
                                 tab_aux(1,1) = (!tau)/100
                                 mtos(tab_aux,ut)
                         endif
                next
next

' So far, we have calculated the p.i.t. series, labelled here as ut , for every model !i and every forecast horizon !h

scalar aux_model = !i
scalar aux_horizon = !h

series zt
for !t = 1 to nobs
               smpl 1985q1+!t-1   1985q1+!t-1
               stom(ut,tab_aux)
               tab_aux = @qnorm(tab_aux)
              mtos(tab_aux,zt)
next

smpl 1985q1 2012q1

equation eq_berkowitz.ls(hac)   zt = c(1) + c(2)*zt(-1)
scalar mu = eq_berkowitz.@coefs(1) / (1-eq_berkowitz.@coefs(2))
scalar rho = eq_berkowitz.@coefs(2)

scalar sig2 = (eq_berkowitz.@se)^2
scalar pi_ = 3.1415926535
scalar T = eq_berkowitz.@regobs 

smpl 1985q1 1985q1
stom(zt,y_aux)
scalar z1 = y_aux(1,1)
scalar sum_aux = 0
matrix(1,1) zt_aux
matrix(1,1) zt1_aux
for !t = 2 to nobs
               smpl 1985q1+!t-2   1985q1+!t-2          'smpl 1984q4+!t-1   1984q4+!t-1   
               stom(zt,zt1_aux)
               smpl 1985q1+!t-1   1985q1+!t-1          'smpl 1984q4+!t       1984q4+!t  
               stom(zt,zt_aux)
               sum_aux = sum_aux + ((  zt_aux(1,1) -mu -rho*zt1_aux(1,1)  )^2) / (2*sig2)
next
scalar log_likelihood1 = -0.5*log(2*pi_) - 0.5*log(sig2/(1-rho^2)) - ((z1-mu/(1-rho))^2) / (2*sig2/(1-rho^2)) - 0.5*(T-1)*log(2*pi_) - 0.5*(T-1)*log(sig2) - sum_aux

if (!h=1) then 
   mu=0    
   rho=0
   sig2=1
   scalar sum_aux = 0
   for !t = 2 to nobs
               smpl 1985q1+!t-2   1985q1+!t-2    '   smpl 1984q4+!t-1   1984q4+!t-1   
               stom(zt,zt1_aux)
               smpl 1985q1+!t-1   1985q1+!t-1     '   smpl 1984q4+!t       1984q4+!t    
               stom(zt,zt_aux)
               sum_aux = sum_aux + ((zt_aux(1,1) -mu -rho*zt1_aux(1,1))^2) / (2*sig2)
   next
   scalar log_likelihood2  = -0.5*log(2*pi_) - 0.5*log(sig2/(1-rho^2)) - ((z1-mu/(1-rho))^2) / (2*sig2/(1-rho^2)) - 0.5*(T-1)*log(2*pi_) - 0.5*(T-1)*log(sig2) - sum_aux
   scalar LR_stat = -2*(log_likelihood2 - log_likelihood1)
   scalar dof = 3   
endif

' Adapted test for longer horizons (i.e., h>1) 
if !h>1 then
   scalar dof = 2    ' Following Jore, Mitchell and Vahey (2008, p.5) and Clements (2004)
   mu=0
   sig2=1
   scalar sum_aux = 0
   for !t = 2 to nobs
               smpl 1985q1+!t-2   1985q1+!t-2    '   smpl 1984q4+!t-1   1984q4+!t-1   
               stom(zt,zt1_aux)
               smpl 1985q1+!t-1   1985q1+!t-1     '   smpl 1984q4+!t       1984q4+!t    
               stom(zt,zt_aux)
               sum_aux = sum_aux + ((zt_aux(1,1) -mu -rho*zt1_aux(1,1))^2) / (2*sig2)
   next
   scalar log_likelihood2  = -0.5*log(2*pi_) - 0.5*log(sig2/(1-rho^2)) - ((z1-mu/(1-rho))^2) / (2*sig2/(1-rho^2)) - 0.5*(T-1)*log(2*pi_) - 0.5*(T-1)*log(sig2) - sum_aux
   scalar LR_stat = -2*(log_likelihood2 - log_likelihood1)
endif

if LR_stat>0 then 
    scalar p_value_Berkowitz = 1-@cchisq(LR_stat,dof) 
endif
if LR_stat<0 or LR_stat=0 then 
    scalar p_value_Berkowitz = 0   
endif
z_p_value_berkowitz(!i,!h) = p_value_Berkowitz


next  'model i
endif  ' h<5, =8, =12
next  'h

endsub
call backtests
