subroutine  backtests

'------------------------------------------------------------------------------------------------------------------------------------------------------------
' ***  Log predictive density scores (LPDS) from Adolfson, Linde and Villani (2005) ***
'The log predictive score is an indicator for the likelihood of several forecasts.
'It is definied as minus the arithmethic mean of the logarithms of the point densities for realized y given data.
matrix(1,1) fan_chart_aux

' Step1: generate for each model !i and horizon !h the conditional quantile series for each out-of-sample period t+h
for !h = h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then

series  y =  pi_aux_h!h                   

for !i = 1 to model_max
series log_pdf_m!i_h!h

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

vector(99) fan_chart

for !tau = 1 to 99
   fan_chart(!tau,1) = 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(!tau) = fan_chart_aux(1,1)

next  ' tau


' Step 2: use a kernel approach to compute the pdf and compare to actual inflation
smpl 1947q1  1971q3   ' 99 observations to convert matrix into an auxiliar  series
mtos(fan_chart,quantiles_y)
matrix(n_grid,2) pdf
matrix(1,1) aux_score

quantiles_y.distdata(dtype=kernel, k=e, ngrid=n_grid) pdf
scalar achou = 0

smpl  1985q1+!t-1   1985q1+!t-1           
   
for !k = 5 to n_grid
          stom(y,y_aux)
          if ( y_aux>pdf(!k,1) ) and (achou=0)  then
                achou = 1       
                aux_score(1,1) = @log10(1e-4)
                if (pdf(!k,2)<>0) then 
                        aux_score(1,1) = @log10(pdf(!k,2))
                endif
                mtos(aux_score,log_pdf_m!i_h!h)
          endif
next   ' k
if (achou=0) then
                aux_score(1,1) = @log10(1e-4)
                mtos(aux_score,log_pdf_m!i_h!h)
endif

next   ' t
next   ' i

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

' end of loopings


' Step 3: compute the average of log scores
matrix(model_max,h_max) z_score
matrix(1,1) aux2_score

for !h=h_min to h_max  
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to model_max
   smpl 1985q1 2012q1      
   series log_aux = log_pdf_m!i_h!h
   z_score(!i,!h) = @mean(log_aux)
next  'model i
endif  ' h<5, =8, =12
next  'h

' end of loopings


' *** Amisano & Giacomini (2007) test ***
matrix(model_max,h_max) z_AG_stat
matrix(model_max,h_max) z_AG_pvalue

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

smpl  1985q1  2012q1               
series benchmark = log_pdf_m35_h!h

for !i = 1 to model_max
    series   diff_lpds_m!i_h!h = benchmark - log_pdf_m!i_h!h
    z_AG_stat(!i,!h) = @mean(diff_lpds_m!i_h!h)
    scalar bw=1.5*!h
   equation Amisano_Giacomini.ls      diff_lpds_m!i_h!h    c
   if (!h=1) then
      if (!i>1) then 
          delete Tab_AG
      endif
      freeze(Tab_AG)   Amisano_Giacomini.ls      diff_lpds_m!i_h!h    c
      z_AG_pvalue(!i,!h) = Tab_AG(9,5)
   endif
   if (!h>1) then
      delete Tab_AG
      freeze(Tab_AG)   Amisano_Giacomini.ls(cov=hac , covbw=bw)  diff_lpds_m!i_h!h    c
      z_AG_pvalue(!i,!h) = Tab_AG(11,5)
   endif

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

next   'h

endsub
call backtests
