'Program written by Bill Clark

close @all

%m=@left(@date,2)					'Current Month, since these are before running, they will reference the real date. In a workfile, it represents an observation date
%d=@mid(@date,4,2)				'Current Day
!year=2000+@right(@date,2)		'Current Year

'------------------------------------------------------------------------------------Estimation Notes -----------------------------------------------------------------------------------

'This program will automatically fit ARIMA models to the data. Note that the max seasonal AR and MA components are NOT the
'same as a seasonal factor. They refer to the seasonal AR and MA components, so they cannot be 12 if monthly. You are NOT
'estimating seasonal factors. Usually an SAR(1) or SMA(1) term is sufficient. You will rarely encounter a genuine SAR(2) or
'SMA(2) process, and even more rarely have enough data to estimate 2 or more seasonal coefficients without the estimation
'algorithm getting into a feedback loop or crashing. Note that if you choose 2 and have monthly data, this means you could have
'2*12=24 seasonal AR terms alone, and then 2*12=24 MA terms if you chose 2 for seasonal MA. This will certainly be a bad
'model due to overfitting and will likley crash most software. EViews will limit TOTAL variables in ARIMA to 24, so choosing
'2 here would not be a good idea. Stick to 1.

'Although a seasonal ARIMA model seems to have only a few parameters, remember
'that backforecasting requires the estimation of one or two seasons' worth of implicit parameters to initialize it.
'Therefore, you should have at least 4 or 5 seasons of data to fit a seasonal ARIMA model. 

'To evaluate a forecast use:

	'@rmse(yf,y)
	'@mae(yf,y)
	'@mape(yf,y)
	'@theil(yf,y)

'ETS model uses automatic model selection


'-------------------------------------------------------------------------- FILL IN THE FOLLOWING --------------------------------------------------------------------------

!annual1=2010  'Historical data start date
!annual2=2018  'This date must be the end of your forecast period

%history_end_date="2017.06"
%forecast_start_date="2017.07"
%forecast_end_date="2018.06"
%wfsaveloc="\\NCRFP4\TAFOA$\EViews\Workfiles\Airports"
%excelsaveloc="\\NCRFP4\TAFOA$\EViews\Excel\Airports\"
!fcst_length=@datediff(@dateval(%forecast_end_date, "yyyy.mm"),@dateval(%forecast_start_date, "yyyy.mm"),"mm")+1

'-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
wfcreate {%wfsaveloc}\tower_ifrvfr_oos_{!year}.{%m}.{%d} u !annual1 !annual2

'Make HWA OOS Page
pagecreate m !annual1 !annual2
pagerename untitled1 hwa_tower_oos
smpl @all

'Get all (*) data and initialize forecast variables
FETCH(d=\\NCRFP4\TAFOA$\EViews\Databases\airport_iv_data_ts_n.edb) *
%towers = "CYVR CYEG CYYC CYXE CYWG CYYZ CYTZ CYOW CYUL CYHZ CYXX CZBB CYBW CYRC CZVL CYFC CYMM CYQX CYHM CYLW CYKF CYNJ CYXU CYQM CYHU CYOO CYPK CYXS CYQB CYQR CYJN CYAM CYYT CYQT CYKZ CYYJ CYXY CYQG CYAV CYZF CYHC"
%flight_rule = "IFR VFR"

table(42,3) a_hwa_feval_ifr
	a_hwa_feval_ifr(1,1) = "Tower"
	a_hwa_feval_ifr(1,2) = "HWA IFR MAPE"
	a_hwa_feval_ifr(1,3) = "HWA IFR MAE"
	a_hwa_feval_ifr.setwidth(2) 14
	a_hwa_feval_ifr.setwidth(3) 14
	a_hwa_feval_ifr.setformat(2,3) f.2

table(42,3) a_hwa_feval_vfr
	a_hwa_feval_vfr(1,1) = "Tower"
	a_hwa_feval_vfr(1,2) = "HWA VFR MAPE"
	a_hwa_feval_vfr(1,3) = "HWA VFR MAE"
	a_hwa_feval_vfr.setwidth(2) 14
	a_hwa_feval_vfr.setwidth(3) 14
	a_hwa_feval_vfr.setformat(2,3) f.2


for %series1 {%flight_rule}
	!row=2
	for %series2 {%towers}
		'Cannot have zeros for any multiplicative analysis, so recode 0 to 0.01
		{%series1}_{%series2}_n = @recode({%series1}_{%series2}_n<=0,0.01,{%series1}_{%series2}_n)
		
		smpl @first %history_end_date
			freeze(mode=overwrite,{%series1}_{%series2}_stats) {%series1}_{%series2}_n.smooth(m) forecast
		smpl %forecast_start_date %forecast_end_date
			genr {%series1}_{%series2}_hwa_f = @round(forecast)

			'Make an actual series that is only for the forecast period
			genr {%series1}_{%series2}_hwa_a = {%series1}_{%series2}_n
			delete forecast
		smpl @all
		
		scalar {%series1}_{%series2}_hwa_mape = @mape({%series1}_{%series2}_hwa_f,{%series1}_{%series2}_hwa_a)
		scalar {%series1}_{%series2}_hwa_mae = @mae({%series1}_{%series2}_hwa_f,{%series1}_{%series2}_hwa_a)

		if %series1="IFR" then		
			a_hwa_feval_ifr(!row,1) = %series2
			a_hwa_feval_ifr(!row,2) = {%series1}_{%series2}_hwa_mape
			a_hwa_feval_ifr(!row,3) = {%series1}_{%series2}_hwa_mae
			!row=!row+1
		else
			a_hwa_feval_vfr(!row,1) = %series2
			a_hwa_feval_vfr(!row,2) = {%series1}_{%series2}_hwa_mape
			a_hwa_feval_vfr(!row,3) = {%series1}_{%series2}_hwa_mae
			!row=!row+1
		endif
	next
next
smpl @all


'Make ETS OOS Page
pagecreate m !annual1 !annual2
pagerename untitled1 ets_tower_oos
smpl @all

'Get all (*) data and initialize forecast variables
FETCH(d=\\NCRFP4\TAFOA$\EViews\Databases\airport_iv_data_ts_n.edb) *
%towers2 = "CYVR CYEG CYYC CYXE CYWG CYYZ CYTZ CYOW CYUL CYHZ CYXX CZBB CYBW CYRC CZVL CYFC CYMM CYQX CYHM CYLW CYKF CYNJ CYXU CYQM CYHU CYOO CYPK CYXS CYQB CYQR CYJN CYAM CYYT CYQT CYKZ CYYJ CYXY CYQG CYAV CYZF CYHC"
%flight_rule2 = "IFR VFR"

table(42,3) a_ets_feval_ifr
	a_ets_feval_ifr(1,1) = "Tower"
	a_ets_feval_ifr(1,2) = "ETS IFR MAPE"
	a_ets_feval_ifr(1,3) = "ETS IFR MAE"
	a_ets_feval_ifr.setwidth(2) 14
	a_ets_feval_ifr.setwidth(3) 14
	a_ets_feval_ifr.setformat(2,3) f.2

table(42,3) a_ets_feval_vfr
	a_ets_feval_vfr(1,1) = "Tower"
	a_ets_feval_vfr(1,2) = "ETS VFR MAPE"
	a_ets_feval_vfr(1,3) = "ETS VFR MAE"
	a_ets_feval_vfr.setwidth(2) 14
	a_ets_feval_vfr.setwidth(3) 14
	a_ets_feval_vfr.setformat(2,3) f.2


for %series1 {%flight_rule2}
	!row=2
	for %series2 {%towers2}
		'Cannot have zeros for any multiplicative analysis, so recode 0 to 0.01
		{%series1}_{%series2}_n = @recode({%series1}_{%series2}_n<=0,0.01,{%series1}_{%series2}_n)
		
		smpl @first %history_end_date
			freeze(mode=overwrite,{%series1}_{%series2}_stats) {%series1}_{%series2}_n.ets(e=e,t=e,s=e,modsel=aic,graph=cl,table=cl) forecast
			
		smpl %forecast_start_date %forecast_end_date
			genr {%series1}_{%series2}_ets_f = @round(forecast)
			
			'Under circumstances of small numbers of observastions (i.e., medium wtc overflights), the ETS framework can make negative forecasts, so recode for this. This is very rare, just check for it
			{%series1}_{%series2}_ets_f  = @recode({%series1}_{%series2}_ets_f <0,0,{%series1}_{%series2}_ets_f )

			'Make an actual series that is only for the forecast period
			genr {%series1}_{%series2}_ets_a = {%series1}_{%series2}_n
			delete forecast
		smpl @all
		
		scalar {%series1}_{%series2}_ets_mape = @mape({%series1}_{%series2}_ets_f,{%series1}_{%series2}_ets_a)
		scalar {%series1}_{%series2}_ets_mae = @mae({%series1}_{%series2}_ets_f,{%series1}_{%series2}_ets_a)

		if %series1="IFR" then		
			a_ets_feval_ifr(!row,1) = %series2
			a_ets_feval_ifr(!row,2) = {%series1}_{%series2}_ets_mape
			a_ets_feval_ifr(!row,3) = {%series1}_{%series2}_ets_mae
			!row=!row+1
		else
			a_ets_feval_vfr(!row,1) = %series2
			a_ets_feval_vfr(!row,2) = {%series1}_{%series2}_ets_mape
			a_ets_feval_vfr(!row,3) = {%series1}_{%series2}_ets_mae
			!row=!row+1
		endif
	next
next
smpl @all


'------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'Clean-up & save
smpl @all
pagedelete untitled
pageselect hwa_tower_oos
wfsave(2) {%wfsaveloc}\tower_ifrvfr_oos_{!year}.{%m}.{%d}.wf1