*
* This is program LIV_RTSIC_roll5_08Nov.prg on 11/21/08
*     It runs real-time forecasts to see if improvement is possible over LIV
*        Compare survey forecast RMSEs to QRTSIC & RTSIC via RMSE
*        Does this with latest available & calculates RMSEs for subsamples
*          Adds subsamples
*       Modified 06Nov to do DM test
*
cal 1947 1 2
alloc 0 2010:2
*
*  set some parameters
*       ffirst = first forecast date
*       flast = last forecast date
*       afirst = first actual date
*       alast = last actual date
*
declare series z
source(noecho) c:\Data\msetest2.src
compute ffirst = 1971:1      ; * first forecast
compute flast = 2007:2       ; * last forecast
compute afirst = 1965:2      ; * first actual
compute alast = 2007:2       ; * last actual
compute aflast = 2006:2      ; * last actual minus one year for lag in forecast
*
open data   c:\data\LIVPIE1.xls
data(format=xls,org=obs) ffirst flast livpie1
*
* First, test capturing the right actuals
*
open data c:\data\LIV_actuals_p_08May.xls
data(format=xls,org=obs,select=2) afirst alast actual_l_q5 actual_1q_q5 actual_1y_q5 $
                                      actual_3y_q5 actual_pb_q5
*
print(dates) afirst flast livpie1 actual_l_q5 actual_1q_q5 actual_1y_q5 $
                                      actual_3y_q5 actual_pb_q5
open data c:\data\SIC_fore_p08.xls
data(format=xls,org=obs,select=2) afirst alast dxq4_hat dxq4q_hat
print(dates) afirst flast livpie1 dxq4_hat dxq4q_hat

compute vtot = 85                                            ;* change this
declare series z1 pi1 pi1a
declare series newfor
declare series pinf
declare vector[series] x(vtot) dxq(vtot) dxq4(vtot)

*  set some parameters
*       send = last observation date in last vintage
*       vend = date of last vintage
*       vtot = number of vintages
*       obstot = number of observations from 1947:2 to send
compute send = 2007:2                                         ;* change this

compute vend = send+1
compute obstot = vtot+75

display 'send = ' send %datelabel(send)
display 'vend = ' vend %datelabel(vend)

display 'vtot = ' vtot
display 'obstot = ' obstot
*
***************************************************
*                                                 *
* Part 5. Calculate RMSEs                         *
*                                                 *
*                                                 *
***************************************************




compute [vector[string]] actualbl = $
    ||'latest','1 quarter','1 year','3 years','pre-bench'||
compute i = 2
*
set z1 = actual_1q_q5
*
display ' '
display 'SIC vs LIV vs SICLA'
display 'Actuals are ' actualbl(i)
display ' '
*
compute flast5 = flast - 3
*
set pi1 / = z1(t+2)
set e    = pi1 - livpie1
set newerror = pi1 - dxq4_hat
set quasie = pi1 - dxq4q_hat
print(dates) ffirst flast5 pi1 livpie1 dxq4_hat e newerror

*
***descriptive statistics of forecast errors

compute naflast = aflast - 11
do fdate = ffirst, naflast
smpl fdate fdate+9
*
*
display ' '
*
***Calculate RMSE
*
compute tse = %nobs
compute aee = %dot(e,e)
compute rmsee = sqrt(aee/tse)
*compute rmse(2,1) = rmsee
display 'This is actuals = ' actualbl(i)
display ' '
display 'Root Mean Square Forecast Error for survey =' rmsee
*
compute tsn = %nobs
compute aen = %dot(newerror,newerror)
compute rmsen = sqrt(aen/tsn)
*compute rmse(2,1) = rmsee
display 'Root Mean Square Forecast Error for SIC forecast =' rmsen
*
set rrmse fdate fdate = rmsen/rmsee
@msetest2(nopictures) livpie1 dxq4_hat pi1 fdate fdate+9 3 p11 p12
*
set prmse fdate fdate = p12


end do fdate
print(dates) ffirst naflast rrmse prmse



*