***forecast correction model for the SPF Survey
*
* This is program SPF_FIEbias_1Qlhs_10yr_psig_08Dec.prg on 12/8/08
*     It runs forecast improvement exercises
*     Uses 1Q actual on LHS of FIE equation
*       Modified to start regressions with at least 20 observations
*       Modified to do DM test
*       Modified to just do FIE when p-value is <.05 for rolling window
*
cal 1947 1 4
alloc 0 2010:4
*
declare series z
declare series z1 pi1 pi1a
declare series newfor
declare series pinf dum
source(noecho) c:\Data\msetest2.src
*
*  set some parameters
*       ffirst = first forecast date
*       flast = last forecast date
*       afirst = first actual date
*       alast = last actual date
*
compute ffirst = 1971:1      ; * first forecast
compute flast = 2008:1       ; * last forecast
compute afirst = 1965:3      ; * first actual
compute alast = 2007:4       ; * last actual
compute aflast = 2006:4      ; * last actual minus one year for lag in forecast

*
open data   c:\data\SPFPIE1_08Jun.xls
data(format=xls,org=obs) ffirst flast spfpie1

open data c:\data\actuals_p_08Mar.xls
data(format=xls,org=obs) afirst alast actual_l_q4 actual_1q_q4 actual_1y_q4 $
                                      actual_3y_q4 actual_pb_q4
*
print(dates) afirst flast spfpie1 actual_l_q4 actual_1q_q4 actual_1y_q4 $
                                      actual_3y_q4 actual_pb_q4
*
*  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:4                                         ;* change this
compute vtot = 170                                            ;* 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
declare vector[series] x(vtot) dxq(vtot) dxq4(vtot)


open data   c:\data\RTDSM_p_08Mar.xls                               ;* change this
data(format=xls,org=obs) 1947:1 send      $
                             P65Q4   $
  P66Q1    P66Q2    P66Q3    P66Q4   $
  P67Q1    P67Q2    P67Q3    P67Q4   $
  P68Q1    P68Q2    P68Q3    P68Q4   $
  P69Q1    P69Q2    P69Q3    P69Q4   $
  P70Q1    P70Q2    P70Q3    P70Q4   $
  P71Q1    P71Q2    P71Q3    P71Q4   $
  P72Q1    P72Q2    P72Q3    P72Q4   $
  P73Q1    P73Q2    P73Q3    P73Q4   $
  P74Q1    P74Q2    P74Q3    P74Q4   $
  P75Q1    P75Q2    P75Q3    P75Q4   $
  P76Q1    P76Q2    P76Q3    P76Q4   $
  P77Q1    P77Q2    P77Q3    P77Q4   $
  P78Q1    P78Q2    P78Q3    P78Q4   $
  P79Q1    P79Q2    P79Q3    P79Q4   $
  P80Q1    P80Q2    P80Q3    P80Q4   $
  P81Q1    P81Q2    P81Q3    P81Q4   $
  P82Q1    P82Q2    P82Q3    P82Q4   $
  P83Q1    P83Q2    P83Q3    P83Q4   $
  P84Q1    P84Q2    P84Q3    P84Q4   $
  P85Q1    P85Q2    P85Q3    P85Q4   $
  P86Q1    P86Q2    P86Q3    P86Q4   $
  P87Q1    P87Q2    P87Q3    P87Q4   $
  P88Q1    P88Q2    P88Q3    P88Q4   $
  P89Q1    P89Q2    P89Q3    P89Q4   $
  P90Q1    P90Q2    P90Q3    P90Q4   $
  P91Q1    P91Q2    P91Q3    P91Q4   $
  P92Q1    P92Q2    P92Q3    P92Q4   $
  P93Q1    P93Q2    P93Q3    P93Q4   $
  P94Q1    P94Q2    P94Q3    P94Q4   $
  P95Q1    P95Q2    P95Q3    P95Q4   $
  P96Q1    P96Q2    P96Q3    P96Q4   $
  P97Q1    P97Q2    P97Q3    P97Q4   $
  P98Q1    P98Q2    P98Q3    P98Q4   $
  P99Q1    P99Q2    P99Q3    P99Q4   $
  P00Q1    P00Q2    P00Q3    P00Q4   $
  P01Q1    P01Q2    P01Q3    P01Q4   $
  P02Q1    P02Q2    P02Q3    P02Q4   $
  P03Q1    P03Q2    P03Q3    P03Q4   $
  P04Q1    P04Q2    P04Q3    P04Q4   $
  P05Q1    P05Q2    P05Q3    P05Q4   $
  P06Q1    P06Q2    P06Q3    P06Q4   $
  P07Q1    P07Q2    P07Q3    P07Q4   $
  P08Q1                                ;* change all vintage names

do v = 1965:4, vend

    compute [label]    vintage1  = %datelabel(v)
    compute [integer]  vintage2  = v - 1965:3

/*
    if v.eq.1965:4
    {
    display @10  'Vintage Index'   @25 'Vintage Date'
    }
    display @10   vintage2         @25 vintage1
*/

end do v




****************************************************
*                                                  *
* Part 1. Index the Vintage and Create             *
*         Annualized Qtr-Over-Qtr Growth Rates     *
*                                                  *
****************************************************
compute [integer] i = 1



dofor z =                                     $
                             P65Q4   $
  P66Q1    P66Q2    P66Q3    P66Q4   $
  P67Q1    P67Q2    P67Q3    P67Q4   $
  P68Q1    P68Q2    P68Q3    P68Q4   $
  P69Q1    P69Q2    P69Q3    P69Q4   $
  P70Q1    P70Q2    P70Q3    P70Q4   $
  P71Q1    P71Q2    P71Q3    P71Q4   $
  P72Q1    P72Q2    P72Q3    P72Q4   $
  P73Q1    P73Q2    P73Q3    P73Q4   $
  P74Q1    P74Q2    P74Q3    P74Q4   $
  P75Q1    P75Q2    P75Q3    P75Q4   $
  P76Q1    P76Q2    P76Q3    P76Q4   $
  P77Q1    P77Q2    P77Q3    P77Q4   $
  P78Q1    P78Q2    P78Q3    P78Q4   $
  P79Q1    P79Q2    P79Q3    P79Q4   $
  P80Q1    P80Q2    P80Q3    P80Q4   $
  P81Q1    P81Q2    P81Q3    P81Q4   $
  P82Q1    P82Q2    P82Q3    P82Q4   $
  P83Q1    P83Q2    P83Q3    P83Q4   $
  P84Q1    P84Q2    P84Q3    P84Q4   $
  P85Q1    P85Q2    P85Q3    P85Q4   $
  P86Q1    P86Q2    P86Q3    P86Q4   $
  P87Q1    P87Q2    P87Q3    P87Q4   $
  P88Q1    P88Q2    P88Q3    P88Q4   $
  P89Q1    P89Q2    P89Q3    P89Q4   $
  P90Q1    P90Q2    P90Q3    P90Q4   $
  P91Q1    P91Q2    P91Q3    P91Q4   $
  P92Q1    P92Q2    P92Q3    P92Q4   $
  P93Q1    P93Q2    P93Q3    P93Q4   $
  P94Q1    P94Q2    P94Q3    P94Q4   $
  P95Q1    P95Q2    P95Q3    P95Q4   $
  P96Q1    P96Q2    P96Q3    P96Q4   $
  P97Q1    P97Q2    P97Q3    P97Q4   $
  P98Q1    P98Q2    P98Q3    P98Q4   $
  P99Q1    P99Q2    P99Q3    P99Q4   $
  P00Q1    P00Q2    P00Q3    P00Q4   $
  P01Q1    P01Q2    P01Q3    P01Q4   $
  P02Q1    P02Q2    P02Q3    P02Q4   $
  P03Q1    P03Q2    P03Q3    P03Q4   $
  P04Q1    P04Q2    P04Q3    P04Q4   $
  P05Q1    P05Q2    P05Q3    P05Q4   $
  P06Q1    P06Q2    P06Q3    P06Q4   $
  P07Q1    P07Q2    P07Q3    P07Q4   $
  P08Q1                                ;* change all vintage names

     set x(i)  /  = z(t)
     compute i = i + 1

end dofor


do i = 1,vtot
   set dxq(i) / = (((x(i)(t)/x(i)(t-1))**4)-1.0)*100.0
end do i






***************************************************
*                                                 *
* Part 2. Replace Missing Values for Various      *
*         Vintages With Data From Preceding       *
*         Vintage                                 *
*                                                 *
*       Everything in this part is                *
*       variable-specific                         *
*                                                 *
***************************************************

* Part 2.A.  Data Sets Feb92 to Nov92
*            Problem:  levels of PGDP exist from 59Q1-present only
*            Solution: replace missing growth rates with those in vintage
*                      Nov91 = x(105)


set dxq(106) 1947:2 1959:1 = dxq(105)(t)
set dxq(107) 1947:2 1959:1 = dxq(105)(t)
set dxq(108) 1947:2 1959:1 = dxq(105)(t)
set dxq(109) 1947:2 1959:1 = dxq(105)(t)

* Part 2.B.  Data Sets Feb96 to Feb97
*            Problem:  levels of PGDP exist from 59Q3-present only
*            Solution: replace missing growth rates with those in vintage
*                      Nov95 = x(121)

set dxq(122) 1947:2 1959:3 = dxq(121)(t)
set dxq(123) 1947:2 1959:3 = dxq(121)(t)
set dxq(124) 1947:2 1959:3 = dxq(121)(t)
set dxq(125) 1947:2 1959:3 = dxq(121)(t)
set dxq(126) 1947:2 1959:3 = dxq(121)(t)

* Part 2.C. Data Set Feb96
*           Problem:  missing 95Q4 observation due to Fed Gov shutdown
*           Solution: replace missing 95Q4 observation with predicted value
*                     from AR(4) specification in growth rates

set x(122) 1995:4 1995:4 = 108.5
set x(122) 1995:3 1995:3 = 107.9
set dxq(122) 1995:4 1995:4 = (((x(122)(t)/x(122)(t-1))**4)-1.0)*100.0
set dxq(122) 1995:3 1995:3 = (((x(122)(t)/x(122)(t-1))**4)-1.0)*100.0
* print 1995:3 1995:4 x(122) dxq(122)

* Part 2.D.  Data Sets Nov99 to Feb00
*            Problem:  levels of PGDP exist from 59Q1-present only
*            Solution: replace missing growth rates with those in vintage
*                      Aug99 = x(136)


set dxq(137) 1947:2 1959:1 = dxq(136)(t)
set dxq(138) 1947:2 1959:1 = dxq(136)(t)



***************************************************
*                                                 *
* Part 3. Construct Indexed Four-Quarter          *
*         Averages                                *
*                                                 *
***************************************************


do i = 1,vtot
*  set dxq4(i) / = ((x(i)(t)/x(i)(t-4))-1.0)*100.0
   set dxq4(i) / = (((((dxq(i)(t)/100.0)+1.0)**0.25)*   $
        (((dxq(i)(t-1)/100.0)+1.0)**0.25)*              $
        (((dxq(i)(t-2)/100.0)+1.0)**0.25)*              $
        (((dxq(i)(t-3)/100.0)+1.0)**0.25))-1.0)*100.0


*  display 'vintage = ' i %datelabel(i+75)
*  print(dates) / dxq(i) dxq4(i)


end do i
*
*
***************************************************
*                                                 *
* Part 4. Create improved forecasts               *
*                                                 *
*                                                 *
***************************************************

*
*
*
*  Here's the new part
*
*  Set up periods for dummy to work
set dum ffirst+28 flast = 0.0
set dum 1989:2 2004:2 = 1.0
*
compute i1 = 50  ; * this is the 1978:1 vintage
*
do fdate = ffirst+28, flast
compute fdate9 = fdate - 9
*
set pinf / = actual_1q_q4(t+4)
* print(dates) ffirst fdate pinf dxq4(i1)
* display 'fdate = ' fdate
* display 'ffirst = ' ffirst
* compute fdate1 = fdate - 1
print(dates) ffirst fdate9 pinf spfpie1

set perror / = pinf - spfpie1

linreg perror ffirst fdate9
# constant
*
set newfor fdate fdate = (dum*%beta(1)) + spfpie1(fdate)
print(dates) fdate fdate newfor spfpie1
display 'alpha=  ' %beta(1)

*
*
compute i1 = i1 + 1
end do fdate
*
***************************************************
*                                                 *
* Part 5. Calculate RMSEs                         *
*                                                 *
*                                                 *
***************************************************





set z1 = actual_1q_q4
*
display ' '
display 'Forecast improvement exercise'
display 'Actuals are 1Q'
display ' '
*
compute ffirst28 = ffirst + 28
compute flast9 = flast - 9
*
set pi1 / = z1(t+4)
set e    = pi1 - spfpie1
set newerror = pi1 - newfor
print(dates) ffirst28 flast9 pi1 spfpie1 newfor e newerror
*
***descriptive statistics of forecast errors

smpl ffirst28 flast9
*
*
display 'Statistics for survey forecast error'
statistics e
*
*
display 'Statistics for improved survey forecast error'
statistics newerror
*
*
display ' '
*
***Calculate RMSE
*
compute ts = %nobs
compute ae = %dot(e,e)
compute rmsee = sqrt(ae/ts)
*compute rmse(2,1) = rmsee
display 'This is actuals = 1Q'
display ' '
display 'Root Mean Square Forecast Error for survey =' rmsee
*
compute ts = %nobs
compute ae = %dot(newerror,newerror)
compute rmsee = sqrt(ae/ts)
*compute rmse(2,1) = rmsee
display 'Root Mean Square Forecast Error for improvement to survey =' rmsee
*
*  End of new part
*
*  Replace this earlier stuff with DM test
*
@msetest2(nopictures) spfpie1 newfor pi1 ffirst28 flast9 5 p11 p1





*