subroutine  backtests

'---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
' RLSW weights (models 1 to 16) from Kascha & Ravazzolo, 2010) based on LPDS: Log predictive density scores, from Adolfson, Linde and Villani (2005) ***
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

smpl @all
series  y =  pi_aux_h!h                   

for !i = 1 to 16
'delete ln_pdf_m!i_h!h
series ln_pdf_m!i_h!h

 for !t = 1 to 165   
 smpl 1970q4+!t   1970q4+!t

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
   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 1970q4+!t   1970q4+!t       
   
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) = @log(1e-4)                
                if (pdf(!k,2)<>0) then 
                        aux_score(1,1) = @log(pdf(!k,2))         
                endif
                mtos(aux_score,ln_pdf_m!i_h!h)
          endif
next   ' k
if (achou=0) then
                aux_score(1,1) = @log(1e-4)                     
                mtos(aux_score,ln_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_weight
matrix(model_max,h_max) z_score
matrix(1,1) aux2_score

for !t = 1 to 109    ' 1985q1 to 2012q1

smpl 1985q1-TS  1985q1+!t-1  ' Evaluate the LPDS weights along the training Sample 
for !h=h_min to h_max
if (!h<5) or (!h=8) or (!h=12) then
for !i = 1 to 16
   series log_aux = ln_pdf_m!i_h!h
   z_score(!i,!h) = @mean(log_aux)
next  'model 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_exp_LPDS = 0
for !i = 1 to 8
    sum_exp_LPDS = sum_exp_LPDS + exp( z_score(!i,!h) )
    series RLSWg1_m!i_h!h 
next  ' !i
for !i = 1 to 8
     z_score_weight(!i,!h) = (exp( z_score(!i,!h) )) / sum_exp_LPDS
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
        aux2_score(1,1) = z_score_weight(!i,!h)
        mtos(aux2_score,RLSWg1_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_exp_LPDS = 0
for !i = 9 to 16
    sum_exp_LPDS = sum_exp_LPDS + exp( z_score(!i,!h) )
    series RLSWg2_m!i_h!h 
next  ' !i
for !i = 9 to 16
     z_score_weight(!i,!h) = (exp( z_score(!i,!h) )) / sum_exp_LPDS
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
        aux2_score(1,1) = z_score_weight(!i,!h)
        mtos(aux2_score,RLSWg2_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_exp_LPDS = 0
for !i = 1 to 16
    sum_exp_LPDS = sum_exp_LPDS + exp( z_score(!i,!h) )
    series RLSWg3_m!i_h!h 
next  ' !i
for !i = 1 to 16
     z_score_weight(!i,!h) = (exp( z_score(!i,!h) )) / sum_exp_LPDS
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
        aux2_score(1,1) = z_score_weight(!i,!h)
        mtos(aux2_score,RLSWg3_m!i_h!h)
next  ' i
endif  ' h<5, =8, =12
next  ' h



next  ' t


endsub
call backtests
