subroutine  qpc
pageselect quarterly

' model 1: AR-GAP (fixed-rho)
' model 2: AR-GAP (2 lags)
' model 3: AR-GAP (RW-AO)
' model 4: PC-GAP (Backward looking)
' model 5: PC-GAP (Backward looking extended)
' model 6: PC-GAP (Forward-Backward looking)
' model 7: PC-GAP (Forward-Backward looking extended)
' model 8: SPF inflation expectation

' Models 1   to 8:   Garch(1,1) Symmetric estimation
' Models 9   to 16: Garch(1,1) Asymmetric estimation (Bank of England)
' Auxiliary models 17 to 24: Quantile regression estimation


' Define below the type of estimation: 
' 1: recursive estimation (always starts in 1961Q1), with expanding sample size      < ----
' 2: recursive estimation (always starts in 1968Q4), with expanding sample size   NOT USED
' 3: rolling estimation(20 year sample)     < ----
' 4: rolling estimation(40 year sample)  NOT USED

scalar estimation = 1   

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

for !h=1 to 12

!TS = 60 ' Training Sample for estimation of combinations of individuals forecasts

pageselect quarterly
smpl @all
scalar TS=!TS
scalar TS_years = @floor(!TS/4 )
scalar TS_quarters = !TS - TS_years*4

if (!h<5) or (!h=8) or (!h=12) then

scalar period_max = 106+!TS + 4           

for !period=1 to period_max    

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

if (estimation=1) then 
   sample sample_estimation             1961q1                                    1984q3+!period+1-!h-!TS
endif
'if (estimation=2) then 
'   sample sample_estimation             1968q4                                    1984q3+!period+1-!h
'endif
if (estimation=3) then 
   sample sample_estimation             1984q3+!period+1-!h-!TS-79      1984q3+!period+1-!h-!TS
endif
'if (estimation=4) then 
'   sample sample_estimation             1984q3+!period+1-!h-159    1984q3+!period+1-!h
'endif

sample sample_forecast_h!h          1984q4+!period-!TS             1984q4+!period-!TS
sample sample_forec_garch_h!h   1984q4+!period+1-!h-!TS    1984q4+!period-!TS
sample sample_last_ten_years      1984q4+!period-40-!TS       1984q4+!period-1-!TS

'----------------------------------------------------------------------------------------------------------------------------------------
pageselect Vintages_GDP_price_index

scalar TS=!TS
scalar TS_years = @floor(!TS/4 )
scalar TS_quarters = !TS - TS_years*4

smpl @all
scalar calc_vint = !period+1-!h +1  ' Following Clark (2011, section 2.Data), since vintage "t" contains information up to "t-1".
' This way, in order to estimate, for instance, the sample 1961q1-1984q4, vintage 85q1 (with information until 84q4) is adopted.
' The only exception is vintage 1996Q1, in which information for GDP price index in 1995Q4 is obtained from next vintage.

if (!h>4) and (!period<20) then       
  scalar calc_vint = !period+1-4+1
endif


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

if (calc_vint=2) then
  scalar year = 1985
  scalar quarter_ = 1
endif

if (calc_vint=3) then
  scalar year = 1985
  scalar quarter_ = 2
endif

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

if (calc_vint>4) then
  scalar calc_years = @floor(calc_vint/4 )
  scalar calc_quarters = calc_vint - calc_years*4
  scalar year = 1984 + calc_years
  scalar quarter_ = 3 + 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

' Adjusting now to training sample
scalar year = year - TS_years
scalar quarter_ = quarter_ - TS_quarters

if (quarter_>4) then
    scalar calc_aux = @floor(quarter_/4)
    scalar quarter_ = quarter_ - calc_aux*4
    scalar year = year + calc_aux
endif
if (quarter_=0) then
    scalar quarter_ = 4
    scalar year = year - 1
endif
if (quarter_=-1) then
    scalar quarter_ = 3
    scalar year = year - 1
endif
if (quarter_=-2) then
    scalar quarter_ = 2
    scalar year = year - 1
endif
if (quarter_=-3) then
    scalar quarter_ = 1
    scalar year = year - 1
endif
if (quarter_=-4) then
    scalar quarter_ = 4
    scalar year = year - 2
endif
if (quarter_=-5) then
    scalar quarter_ = 3
    scalar year = year - 2
endif
if (quarter_=-6) then
    scalar quarter_ = 2
    scalar year = year - 2
endif
if (quarter_=-7) then
    scalar quarter_ = 1
    scalar year = year - 2
endif
if (quarter_=-8) then
    scalar quarter_ = 4
    scalar year = year - 3
endif
if (quarter_=-9) then
    scalar quarter_ = 3
    scalar year = year - 3
endif
if (quarter_=-10) then
    scalar quarter_ = 2
    scalar year = year - 3
endif
if (quarter_=-11) then
    scalar quarter_ = 1
    scalar year = year - 3
endif
if (quarter_=-12) then
    scalar quarter_ = 4
    scalar year = year - 4
endif

for !i=year to year
for !q= quarter_ to quarter_
   series gdp_price_vintage = gdp_price_p!iq!q 
next
next

' No real time data with the two lines below activated
'smpl @all
'series gdp_price_vintage = gdp_price_p2011q2

'--------------------------------------------------------------------------------------------------------
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))



delete inf_trend
'smpl     1984q4+!period-40-!TS       1984q4+!period-1-!TS  ' last ten  years = 40 observ.
'smpl     1984q4+!period-20-!TS       1984q4+!period-1-!TS  ' last five years = 20 observ.
'equation eq_inf_trend.ls    pi_a   mean_cpi10
'smpl     1947q1   1984q4+!period-1-!TS
'series inf_trend_aux = mean_cpi10 * eq_inf_trend.@coefs(1)
smpl   1947q1     1966q4
series inf_trend = mean_cpi10
smpl  1967q1      1984q3+!period+1-!h-!TS
series inf_trend = @movav(pi_a,16)       ' @movav(pi_a,20)
smpl     1984q3+!period+1-!h-!TS+1     2016q4   ' forecast of inflation trend = random walk
series inf_trend = inf_trend(-1)


smpl @all
if (!h<4) then
   series pi = pi_a
   series gap = pi - inf_trend
endif
if (!h=4) then
   series pi = pi_b
   series gap = pi - inf_trend
endif
if (!h>4) then
   series pi = pi_b
   series gap = pi - inf_trend
endif

'--------------------------------------------------------------------------------------------------------
' Point (and density) forecast definition is f ( t ; t-h)  
' Expectation is defined as exp_inf_t+h ; t. 
' Since y(t,h) = log(1+gap(t)/100) for all horizon (h), then, expectation as covariate is used as exp_inf_t+h;t (-h)=exp_inf_t;t-h (0), i.e., lagged h periods.

series y_h1   =   log(1+gap/100)
series y_h2   =   log(1+gap/100)
series y_h3   =   log(1+gap/100)
series y_h4   =   log(1+gap/100)
series y_h8   =   log(1+gap/100)
series y_h12   =   log(1+gap/100)

group  x1_h1   log(1+gap(-1)/100)     
group  x1_h2   log(1+gap(-2)/100)    
group  x1_h3   log(1+gap(-3)/100)     
group  x1_h4   log(1+gap(-4)/100)       
group  x1_h8   log(1+gap(-8)/100)       
group  x1_h12   log(1+gap(-12)/100)       

group  x2_h1    one      log(1+gap(-1)/100)      log(1+gap(-2)/100)   
group  x2_h2    one      log(1+gap(-2)/100)      log(1+gap(-3)/100)   
group  x2_h3    one      log(1+gap(-3)/100)      log(1+gap(-4)/100)   
group  x2_h4    one      log(1+gap(-4)/100)      log(1+gap(-5)/100)   
group  x2_h8    one      log(1+gap(-8)/100)      log(1+gap(-9)/100)   
group  x2_h12  one      log(1+gap(-12)/100)   log(1+gap(-13)/100)   

group  x3_h1     (log(1+gap(-1)/100) + log(1+gap(-2)/100) + log(1+gap(-3)/100) + log(1+gap(-4)/100))     
group  x3_h2     (log(1+gap(-2)/100) + log(1+gap(-3)/100) + log(1+gap(-4)/100) + log(1+gap(-5)/100))     
group  x3_h3     (log(1+gap(-3)/100) + log(1+gap(-4)/100) + log(1+gap(-5)/100) + log(1+gap(-6)/100))     
group  x3_h4     (log(1+gap(-4)/100) + log(1+gap(-5)/100) + log(1+gap(-6)/100) + log(1+gap(-7)/100))     
group  x3_h8     (log(1+gap(-8)/100) + log(1+gap(-9)/100) + log(1+gap(-10)/100) + log(1+gap(-11)/100))     
group  x3_h12   (log(1+gap(-12)/100) + log(1+gap(-13)/100) + log(1+gap(-14)/100) + log(1+gap(-15)/100))     

group  x4_h1     one      log(1+gap(-1)/100)        log(1+unemp(-1)/100)   
group  x4_h2     one      log(1+gap(-2)/100)        log(1+unemp(-2)/100)   
group  x4_h3     one      log(1+gap(-3)/100)        log(1+unemp(-3)/100)   
group  x4_h4     one      log(1+gap(-4)/100)        log(1+unemp(-4)/100)   
group  x4_h8     one      log(1+gap(-8)/100)        log(1+unemp(-8)/100)   
group  x4_h12   one      log(1+gap(-12)/100)      log(1+unemp(-12)/100)   

group  x5_h1     one      log(1+gap(-1)/100)       log(1+gap(-2)/100)      log(1+unemp(-1)/100)     log(1+unemp(-2)/100)   
group  x5_h2     one      log(1+gap(-2)/100)       log(1+gap(-3)/100)      log(1+unemp(-2)/100)     log(1+unemp(-3)/100)   
group  x5_h3     one      log(1+gap(-3)/100)       log(1+gap(-4)/100)      log(1+unemp(-3)/100)     log(1+unemp(-4)/100)   
group  x5_h4     one      log(1+gap(-4)/100)       log(1+gap(-5)/100)      log(1+unemp(-4)/100)     log(1+unemp(-5)/100)   
group  x5_h8     one      log(1+gap(-8)/100)       log(1+gap(-9)/100)      log(1+unemp(-8)/100)     log(1+unemp(-9)/100)   
group  x5_h12   one      log(1+gap(-12)/100)       log(1+gap(-13)/100)      log(1+unemp(-12)/100)     log(1+unemp(-13)/100)   

group  x6_h1     one      log(1+(exp_gdp_inf_h1(-1)-inf_trend(-1))/100)             log(1+gap(-1)/100)      log(1+unemp(-1)/100)    
group  x6_h2     one      log(1+(exp_gdp_inf_h2(-2)-inf_trend(-2))/100)             log(1+gap(-2)/100)      log(1+unemp(-2)/100)    
group  x6_h3     one      log(1+(exp_gdp_inf_h3(-3)-inf_trend(-3))/100)             log(1+gap(-3)/100)      log(1+unemp(-3)/100)    
group  x6_h4     one      log(1+(exp_gdp_inf_h4(-4)-inf_trend(-4))/100)             log(1+gap(-4)/100)      log(1+unemp(-4)/100)    
group  x6_h8     one      log(1+(exp_gdp_inf_h8(-8)-inf_trend(-8))/100)             log(1+gap(-8)/100)      log(1+unemp(-8)/100)    
group  x6_h12   one      log(1+(exp_gdp_inf_h12(-12)-inf_trend(-12))/100)      log(1+gap(-12)/100)   log(1+unemp(-12)/100)    

group  x7_h1     one      log(1+(exp_gdp_inf_h1(-1)-inf_trend(-1))/100)             log(1+gap(-1)/100)      log(1+unemp(-1)/100)     dlog(1+unemp(-1)/100)    
group  x7_h2     one      log(1+(exp_gdp_inf_h2(-2)-inf_trend(-2))/100)             log(1+gap(-2)/100)      log(1+unemp(-2)/100)     dlog(1+unemp(-2)/100)    
group  x7_h3     one      log(1+(exp_gdp_inf_h3(-3)-inf_trend(-3))/100)             log(1+gap(-3)/100)      log(1+unemp(-3)/100)     dlog(1+unemp(-3)/100)    
group  x7_h4     one      log(1+(exp_gdp_inf_h4(-4)-inf_trend(-4))/100)             log(1+gap(-4)/100)      log(1+unemp(-4)/100)     dlog(1+unemp(-4)/100)    
group  x7_h8     one      log(1+(exp_gdp_inf_h8(-8)-inf_trend(-8))/100)             log(1+gap(-8)/100)      log(1+unemp(-8)/100)     dlog(1+unemp(-8)/100)    
group  x7_h12   one      log(1+(exp_gdp_inf_h12(-12)-inf_trend(-12))/100)      log(1+gap(-12)/100)    log(1+unemp(-12)/100)  dlog(1+unemp(-12)/100)    

'**************************************************************************
' Part I - Initialization 

scalar mu
scalar sig
scalar sig2
scalar skw
scalar skwn
scalar a_aux
scalar s_aux
scalar s1
scalar s2
scalar s3
scalar s4
scalar aux1
scalar aux2
scalar sig_1
scalar sig_2
scalar C_aux
scalar p
scalar pi_ = 3.1415926535

matrix(1,1) mu_aux
matrix(1,1) sig_aux
matrix(1,1) skw_aux
matrix(1,1) aux_garch

for !m=1 to 16
series model!m_mu  
series model!m_mu_  
series model!m_se
series model!m_se_
series model!m_sig2
series model!m_sig2_
series model!m_skw
series model!m_skw_norm
series model!m_check
next

'--------------------------------------------------------------------------------------------------
' Garch estimation (initialization part):  Models 1 to 8
'--------------------------------------------------------------------------------------------------

' --------------------------------------- Model 1 ----------------------------------------------
scalar rho = 0.46
scalar rho_h!h = rho^!h

smpl  sample_estimation
equation model1_h!h.arch(1,1)  y_h!h   =    rho_h!h *  log(1+gap(-!h)/100)     

if (!period>90) then
scalar aux_ncoef = model1_h!h.@ncoef
scalar aux_coef1 = model1_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model1_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model1_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then  
   equation model1_h!h.arch(1,1,vt)  y_h!h   =    rho_h!h *  log(1+gap(-!h)/100)     
endif
endif

smpl  sample_forec_garch_h!h
model1_h!h.forecast(f=na)    model1_mu_  model1_se_   model1_sig2_ 
smpl  sample_forecast_h!h
model1_mu  = model1_mu_  
model1_se = 1e-6
stom(model1_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model1_se   = (model1_sig2_)^0.5
endif
' --------------------------------------- Model 2 ----------------------------------------------
smpl  sample_estimation
equation model2_h!h.arch(1,1)  y_h!h   x2_h!h

if (!period>90) then
scalar aux_ncoef = model2_h!h.@ncoef
scalar aux_coef1 = model2_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model2_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model2_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then  
   equation model2_h!h.arch(1,1,vt)   y_h!h   x2_h!h 
endif
endif

smpl  sample_forec_garch_h!h
model2_h!h.forecast(f=na)    model2_mu_  model2_se_   model2_sig2_ 
smpl  sample_forecast_h!h
model2_mu  = model2_mu_  
model2_se = 1e-6
stom(model2_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model2_se   = (model2_sig2_)^0.5
endif
' --------------------------------------- Model 3 ----------------------------------------------
smpl  sample_estimation
equation model3_h!h.arch(1,1)  y_h!h   =   0.25 *  (log(1+gap(-!h)/100) + log(1+gap(-!h-1)/100) + log(1+gap(-!h-2)/100) + log(1+gap(-!h-3)/100))    

if (!period>90) then
scalar aux_ncoef = model3_h!h.@ncoef
scalar aux_coef1 = model3_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model3_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model3_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then  
   equation model3_h!h.arch(1,1,vt)   y_h!h   =   0.25 *  (log(1+gap(-!h)/100) + log(1+gap(-!h-1)/100) + log(1+gap(-!h-2)/100) + log(1+gap(-!h-3)/100))  
endif
endif

smpl  sample_forec_garch_h!h
model3_h!h.forecast(f=na)    model3_mu_  model3_se_   model3_sig2_ 
smpl  sample_forecast_h!h
model3_mu  = model3_mu_  
model3_se = 1e-6
stom(model3_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model3_se   = (model3_sig2_)^0.5
endif
' --------------------------------------- Model 4 ----------------------------------------------
smpl  sample_estimation
equation model4_h!h.arch(1,1)  y_h!h   x4_h!h

if (!period>90) then
scalar aux_ncoef = model4_h!h.@ncoef
scalar aux_coef1 = model4_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model4_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model4_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then   
   equation model4_h!h.arch(1,1,vt)   y_h!h   x4_h!h 
endif
endif

smpl  sample_forec_garch_h!h
model4_h!h.forecast(f=na)    model4_mu_  model4_se_   model4_sig2_ 
smpl  sample_forecast_h!h
model4_mu  = model4_mu_  
model4_se = 1e-6
stom(model4_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model4_se   = (model4_sig2_)^0.5
endif
' --------------------------------------- Model 5 ----------------------------------------------
smpl  sample_estimation
equation model5_h!h.arch(1,1)  y_h!h   x5_h!h

if (!period>90) then
scalar aux_ncoef = model5_h!h.@ncoef
scalar aux_coef1 = model5_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model5_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model5_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then   
   equation model5_h!h.arch(1,1,vt)   y_h!h   x5_h!h 
endif
endif

smpl  sample_forec_garch_h!h
model5_h!h.forecast(f=na)    model5_mu_  model5_se_   model5_sig2_ 
smpl  sample_forecast_h!h
model5_mu  = model5_mu_  
model5_se = 1e-6
stom(model5_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model5_se   = (model5_sig2_)^0.5
endif
' --------------------------------------- Model 6 ----------------------------------------------
smpl  sample_estimation
equation model6_h!h.arch(1,1)  y_h!h   x6_h!h

if (!period>90) then
scalar aux_ncoef = model6_h!h.@ncoef
scalar aux_coef1 = model6_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model6_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model6_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then   
   equation model6_h!h.arch(1,1,vt)   y_h!h   x6_h!h 
endif
endif

smpl  sample_forec_garch_h!h
model6_h!h.forecast(f=na)    model6_mu_  model6_se_   model6_sig2_ 
smpl  sample_forecast_h!h
model6_mu  = model6_mu_  
model6_se = 1e-6
stom(model6_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model6_se   = (model6_sig2_)^0.5
endif
' --------------------------------------- Model 7 ----------------------------------------------
smpl  sample_estimation
equation model7_h!h.arch(1,1)  y_h!h   x7_h!h

if (!period>90) then
scalar aux_ncoef = model7_h!h.@ncoef
scalar aux_coef1 = model7_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model7_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model7_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then   
   equation model7_h!h.arch(1,1,vt)   y_h!h   x7_h!h 
endif
endif

smpl  sample_forec_garch_h!h
model7_h!h.forecast(f=na)    model7_mu_  model7_se_   model7_sig2_ 
smpl  sample_forecast_h!h
model7_mu  = model7_mu_  
model7_se = 1e-6
stom(model7_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model7_se   = (model7_sig2_)^0.5
endif
' --------------------------------------- Model 8 ----------------------------------------------
smpl  sample_estimation
equation model8_h!h.arch(1,1)  y_h!h    log(1+(exp_gdp_inf_h!h(-!h)-inf_trend(-!h))/100)

if (!period>90) then
scalar aux_ncoef = model8_h!h.@ncoef
scalar aux_coef1 = model8_h!h.@coefs(aux_ncoef-2)
scalar aux_coef2 = model8_h!h.@coefs(aux_ncoef-1)
scalar aux_coef3 = model8_h!h.@coefs(aux_ncoef)
if  ((aux_coef2+aux_coef3)>1) or (aux_coef1<0) then   
   equation model8_h!h.arch(1,1,vt)   y_h!h    log(1+(exp_gdp_inf_h!h(-!h)-inf_trend(-!h))/100) 
endif
endif

smpl  sample_forec_garch_h!h
model8_h!h.forecast(f=na)    model8_mu_  model8_se_   model8_sig2_ 
smpl  sample_forecast_h!h
model8_mu  = model8_mu_  
model8_se = 1e-6
stom(model8_sig2_,aux_garch)
if (aux_garch(1,1)>0) then 
   model8_se   = (model8_sig2_)^0.5
endif

'**********************************************************
' Part II - Density forecast (out-of-sample) 
'**********************************************************

matrix(9,1) coefs9
matrix(7,1) coefs7
matrix(6,1) coefs6
matrix(5,1) coefs5
matrix(4,1) coefs4
matrix(3,1) coefs3
matrix(2,1) coefs2

matrix(1,9) x_aux9
matrix(1,7) x_aux7
matrix(1,6) x_aux6
matrix(1,5) x_aux5
matrix(1,4) x_aux4
matrix(1,3) x_aux3
matrix(1,2) x_aux2
matrix(1,1) x_aux1
matrix(1,1) forecast

for !tau=1 to 99

'------------------------------------------------------------------------------------------------------
' Model 1 
series forec
smpl  sample_forecast_h!h
forec =  model1_mu + model1_se*@qnorm(!tau/100)
series m1_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 2 
series forec
smpl  sample_forecast_h!h
forec =  model2_mu + model2_se*@qnorm(!tau/100)
series m2_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 3 
series forec
smpl  sample_forecast_h!h
forec =  model3_mu + model3_se*@qnorm(!tau/100)
series m3_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 4 
series forec
smpl  sample_forecast_h!h
forec =  model4_mu + model4_se*@qnorm(!tau/100)
series m4_h!h_q!tau =  ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 5 
series forec
smpl  sample_forecast_h!h
forec =  model5_mu + model5_se*@qnorm(!tau/100)
series m5_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 6 
series forec
smpl  sample_forecast_h!h
forec =  model6_mu + model6_se*@qnorm(!tau/100)
series m6_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 7 
series forec
smpl  sample_forecast_h!h
forec =  model7_mu + model7_se*@qnorm(!tau/100)
series m7_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 8 
series forec
smpl  sample_forecast_h!h
forec =  model8_mu + model8_se*@qnorm(!tau/100)
series m8_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 



'------------------------------------------------------------------------------------------------------
' Formulas below from John (1982) "The three-parameter two-piece normal family of distributions and its fitting" and Julio (2006) "The fan chart: The technical details of the new implementation".
'------------------------------------------------------------------------------------------------------
' Model 9  (equal to model 1 but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model1_mu,mu_aux,sample_forecast_h!h)
stom(model1_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model9_mu = mu
model9_sig2 = sig2
model9_skw = skw
model9_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m9_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 10  (equal to model 2 but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model2_mu,mu_aux,sample_forecast_h!h)
stom(model2_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model10_mu = mu
model10_sig2 = sig2
model10_skw = skw
model10_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m10_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 11  (equal to model 3 but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model3_mu,mu_aux,sample_forecast_h!h)
stom(model3_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model11_mu = mu
model11_sig2 = sig2
model11_skw = skw
model11_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m11_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 12  (equal to model 4 but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model4_mu,mu_aux,sample_forecast_h!h)
stom(model4_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model12_mu = mu
model12_sig2 = sig2
model12_skw = skw
model12_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m12_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 13 (equal to model 5  but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model5_mu,mu_aux,sample_forecast_h!h)
stom(model5_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model13_mu = mu
model13_sig2 = sig2
model13_skw = skw
model13_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m13_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 14 (equal to model 6  but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model6_mu,mu_aux,sample_forecast_h!h)
stom(model6_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model14_mu = mu
model14_sig2 = sig2
model14_skw = skw
model14_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m14_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 15 (equal to model 7  but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model7_mu,mu_aux,sample_forecast_h!h)
stom(model7_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model15_mu = mu
model15_sig2 = sig2
model15_skw = skw
model15_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m15_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------
' Model 16 (equal to model 8  but with asymmetric Garch density estimation)   
smpl sample_last_ten_years
skw = @skew(gap)   ' If skw is set to zero then model becomes a standard Garch
skwn = skw / (1+((skw^2)^0.5))
smpl  sample_forecast_h!h

stom(model8_mu,mu_aux,sample_forecast_h!h)
stom(model8_se,sig_aux,sample_forecast_h!h)
mu = mu_aux(1,1)
sig = sig_aux(1,1)
sig2 = sig^2
model16_mu = mu
model16_sig2 = sig2
model16_skw = skw
model16_skw_norm = skwn

sig_1 = (sig2/(1-skwn))^0.5
sig_2 = (sig2/(1+skwn))^0.5
C_aux = ((2/pi_)^0.5) / (sig_1+sig_2)
p = sig_1 / (sig_1 + sig_2)
if (!tau/100<p) then
forecast = mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100=p) then
forecast =  mu + sig_1 * @qnorm(   (!tau/100) / (C_aux*sig_1*((2*pi_)^0.5))    )
endif
if (!tau/100>p) then
forecast =  mu + sig_2 * @qnorm(   (!tau/100 - 1 + C_aux*sig_2*((2*pi_)^0.5)) / (C_aux*sig_2*((2*pi_)^0.5))    )
endif
series forec
mtos(forecast,forec)
series m16_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

'------------------------------------------------------------------------------------------------------
' Models 17-24: QR estimation
'------------------------------------------------------------------------------------------------------
' Model 17  (equal to model 1 but with QR estimation)   
scalar rho = 0.46
scalar rho_h!h = rho^!h
smpl @all
series temp17 =  y_h!h - rho_h!h * log(1+gap(-!h)/100)     
smpl  sample_estimation
equation model17_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)  temp17  c
smpl  sample_forecast_h!h
series temp17b   = rho_h!h * log(1+gap(-!h)/100)     
stom(temp17b,aux_c1)
forecast(1,1) = model17_h!h.@coef(1)  +  aux_c1(1,1)
series forec
mtos(forecast,forec)
series m17_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 18  (equal to model 2 but with QR estimation)   
smpl  sample_estimation
equation model18_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)   y_h!h     c      log(1+gap(-!h)/100)        log(1+gap(-!h-1)/100)     
smpl  sample_forecast_h!h
series temp2   = log(1+gap(-!h)/100)
series temp3   = log(1+gap(-!h-1)/100)
stom(temp2,aux_c2)
stom(temp3,aux_c3)
forecast(1,1) = model18_h!h.@coef(1)  +  model18_h!h.@coef(2) * aux_c2(1,1)  +  model18_h!h.@coef(3) * aux_c3(1,1)
series forec
mtos(forecast,forec)
series m18_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 19  (equal to model 3 but with QR estimation)   
smpl @all
series temp19a = 0.25 * (log(1+gap(-!h)/100) + log(1+gap(-!h-1)/100) + log(1+gap(-!h-2)/100) + log(1+gap(-!h-3)/100))  
series temp19b = y_h!h  -  temp19a
smpl  sample_estimation
equation model19_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)   temp19b  c     
smpl  sample_forecast_h!h
stom(temp19a,aux_c1)
forecast(1,1) = model19_h!h.@coef(1) + aux_c1(1,1)
series forec
mtos(forecast,forec)
series m19_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 20  (equal to model 4 but with QR estimation)   
smpl  sample_estimation
equation model20_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)   y_h!h     c      log(1+gap(-!h)/100)        log(1+unemp(-!h)/100)   
smpl  sample_forecast_h!h
series temp2   = log(1+gap(-!h)/100)
series temp3   = log(1+unemp(-!h)/100)   
stom(temp2,aux_c2)
stom(temp3,aux_c3)
forecast(1,1) = model20_h!h.@coef(1)  +  model20_h!h.@coef(2) * aux_c2(1,1)  +  model20_h!h.@coef(3) * aux_c3(1,1)
series forec
mtos(forecast,forec)
series m20_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 21  (equal to model 5 but with QR estimation)   
smpl  sample_estimation
equation model21_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)   y_h!h     c      log(1+gap(-!h)/100)     log(1+gap(-!h-1)/100)       log(1+unemp(-!h)/100)     log(1+unemp(-!h-1)/100)   
smpl  sample_forecast_h!h
series temp2   = log(1+gap(-!h)/100)
series temp3   = log(1+gap(-!h-1)/100)
series temp4   = log(1+unemp(-!h)/100)   
series temp5   = log(1+unemp(-!h-1)/100)   
stom(temp2,aux_c2)
stom(temp3,aux_c3)
stom(temp4,aux_c4)
stom(temp5,aux_c5)
forecast(1,1) = model21_h!h.@coef(1)  +  model21_h!h.@coef(2) * aux_c2(1,1)  +  model21_h!h.@coef(3) * aux_c3(1,1)  +  model21_h!h.@coef(4) * aux_c4(1,1)  +  model21_h!h.@coef(5) * aux_c5(1,1)
series forec
mtos(forecast,forec)
series m21_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 22  (equal to model 6 but with QR estimation)   
smpl @all
series expec = log(1+(exp_gdp_inf_h{!h}(-!h)-inf_trend(-!h))/100)
smpl  sample_estimation
equation model22_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)   y_h!h     c     expec     log(1+gap(-!h)/100)     log(1+unemp(-!h)/100)    
smpl  sample_forecast_h!h
series temp2   = expec
series temp3   = log(1+gap(-!h)/100)
series temp4   = log(1+unemp(-!h)/100)   
stom(temp2,aux_c2)
stom(temp3,aux_c3)
stom(temp4,aux_c4)
forecast(1,1) = model22_h!h.@coef(1)  +  model22_h!h.@coef(2) * aux_c2(1,1)  +  model22_h!h.@coef(3) * aux_c3(1,1)  +  model22_h!h.@coef(4) * aux_c4(1,1)  
series forec
mtos(forecast,forec)
series m22_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 23  (equal to model 7 but with QR estimation)   
smpl @all
series expec = log(1+(exp_gdp_inf_h{!h}(-!h)-inf_trend(-!h))/100)
smpl  sample_estimation
equation model23_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)   y_h!h     c     expec     log(1+gap(-!h)/100)     log(1+unemp(-!h)/100)      dlog(1+unemp(-!h)/100)    
smpl  sample_forecast_h!h
series temp2   = expec
series temp3   = log(1+gap(-!h)/100)
series temp4   = log(1+unemp(-!h)/100)   
series temp5   = dlog(1+unemp(-!h)/100)   
stom(temp2,aux_c2)
stom(temp3,aux_c3)
stom(temp4,aux_c4)
stom(temp5,aux_c5)
forecast(1,1) = model23_h!h.@coef(1)  +  model23_h!h.@coef(2) * aux_c2(1,1)  +  model23_h!h.@coef(3) * aux_c3(1,1)  +  model23_h!h.@coef(4) * aux_c4(1,1)  +  model23_h!h.@coef(5) * aux_c5(1,1)  
series forec
mtos(forecast,forec)
series m23_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

' Model 24  (equal to model 8 but with QR estimation)   
smpl @all
series expec = log(1+(exp_gdp_inf_h{!h}(-!h)-inf_trend(-!h))/100)
series temp24 =  y_h!h   -   expec   
smpl  sample_estimation
equation model24_h!h.qreg(quant =  !tau/100 ,  cov=sandwich)    temp24  c
smpl  sample_forecast_h!h
series temp   = expec
stom(temp,aux_c1)
forecast(1,1) = model24_h!h.@coef(1) + aux_c1(1,1) 
series forec
mtos(forecast,forec)
series m24_h!h_q!tau =   ((exp( forec ) - 1 ) * 100) + inf_trend(-!h) 
'------------------------------------------------------------------------------------------------------

next  ' !tau

delete sample_estimation
delete sample_forecast_h!h
delete sample_forec_garch_h!h  
delete  sample_last_ten_years
'delete  sample_gap

next  ' !period

endif 'h<5, =8 or =12
next  ' forecast horizon (h)

smpl @all

endsub
call qpc
