' THIS PROGRAM CALCULATES DYNAMIC FACTOR MODELS AND USES THEM TO FORECAST.  THE DFMs ARE CALCULATED USING VARIOUS METHODS THAT ARE SUMMARIZED IN THE SECTION.

%prog_dir = "C:\Arora_Files\EIA_Files\Projects\Other\2013\Factor_Models\Data"
'Set main directory

%data_page = "ngdata"
' Page with the data

%dfm_page = "dfm_g2"
' Page to work on.

%reg = "ak gm ls nm ok os tx wy"
%all = "ak gm ls nm ok os tx wy us"
%reg_apc = "ak_apc gm_apc ls_apc nm_apc ok_apc os_apc tx_apc wy_apc"
%all_apc = "ak_apc gm_apc ls_apc nm_apc ok_apc os_apc tx_apc wy_apc us_apc"
' Gather production regions for use in various loops

!max_lag = 9
' Specify the number of lags to use of independent variables in equations: max is 9.

'pageselect {%data_page}
'!data_begin = @dateval(@otod(1))
'!data_end = @dateval("201305","YYYYMM")
'' Set first and last data points

!exist = @pageexist(%dfm_page)
	if !exist = 1 then
	pagedelete {%dfm_page}
	pagecreate(page={%dfm_page}) m 1973m01 2013m05
	else
	pagecreate(page={%dfm_page}) m 1973m01 2013m05
	endif
'Create page or delete then create

'%%% BEGIN WORKING ON MAIN PAGE

pageselect {%dfm_page}
' Select page with data

for %y {%all}
	copy {%data_page}\{%y}_ngpd {%dfm_page}\{%y}
next
' Copy over data series and rename

group g_reg_raw {%reg}
group g_all_raw {%all}
' Gather all regions and total us production with raw data

for %y {%all}
	formula {%y}_apc = @pcy({%y})
'	{%y}_apc.trim(0.10,0.90,wins) _trim
'	delete {%y}_apc
'	genr {%y}_apc = {%y}_apc_trim
'	delete {%y}_apc_trim
next
'' Get annual percent changes, but also get rid of major outliers using winsoring.

group g_reg_apc {%reg_apc}
group g_all_apc {%all_apc}
' Gather all regions and total us production with raw data

'%%% BEGIN DFM ANALYSIS USING METHOD OF SW (2002)

smpl 2007m01 @last
' Only use data after shale

for %y reg all
	g_{%y}_apc.pcomp(eigval=reg_ev_cov,cov=cov)
	freeze({%y}_sc_cov) g_{%y}_apc.pcomp(cov=cov,out=graph,scree)

	g_{%y}_apc.pcomp(eigval=reg_ev_cor,cov=corr)
	freeze({%y}_sc_cor) g_{%y}_apc.pcomp(cov=corr,out=graph,scree)
next
' Write out eigenvalues and save scree plot for each group using the covariance or correlation matrix

for %y reg all
	for %t cov cor
		for !j = 2 to 5
		
		if %t = "cov" then
			factor f_!j_{%y}_{%t}.pf(n=!j,priors=frac,priorfrac=1,cov=cov) g_{%y}_apc
		else
		if %t = "cor" then
			factor f_!j_{%y}_{%t}.pf(n=!j,priors=frac,priorfrac=1,cov=corr) g_{%y}_apc
		endif
		endif
		' Estimate factor model using PCA (priorfrac=1) and working on either covariance or correlation matrix.

		f_!j_{%y}_{%t}.rotate(type=orthog, method=varimax)
		' Rotate factors.
		
		if !j = 2 then
			f_!j_{%y}_{%t}.makescores(n=fs_!j_{%y}_{%t}) fs_!j1_{%y}_{%t} fs_!j2_{%y}_{%t}
		else
		if !j = 3 then
			f_!j_{%y}_{%t}.makescores(n=fs_!j_{%y}_{%t}) fs_!j1_{%y}_{%t} fs_!j2_{%y}_{%t} fs_!j3_{%y}_{%t}
		else
		if !j = 4 then
			f_!j_{%y}_{%t}.makescores(n=fs_!j_{%y}_{%t}) fs_!j1_{%y}_{%t} fs_!j2_{%y}_{%t} fs_!j3_{%y}_{%t} fs_!j4_{%y}_{%t}
		else
		if !j = 5 then
			f_!j_{%y}_{%t}.makescores(n=fs_!j_{%y}_{%t}) fs_!j1_{%y}_{%t} fs_!j2_{%y}_{%t} fs_!j3_{%y}_{%t} fs_!j4_{%y}_{%t} fs_!j5_{%y}_{%t}
		endif
		endif
		endif
		endif
		' Back out scores
		next		
	next
next
' Generate the factor model for each group using the covariance or correlation matrix and then rotate for clearer interpretation.
' Also generate scores using rotated loadings and standard regression and put into separate series.  Even though the factors do not show it, these work on rotated scores.
' Notice that we also loop over the number of factors, mainly for sensitivity analysis.

close @objects
' Close any open objects

'%%% SET UP AND ESTIMATE REGRESSIONS USING FACTORS

for %y reg all
	for %t cov cor
		for !j = 2 to 5
			for !k = 1 to !max_lag
			!k6 = -5-!k
			!k12 = -11-!k
			if !j = 2 then
				equation eq!ka_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-1)) fs_!j1_{%y}_{%t}(-1 to -!k) fs_!j2_{%y}_{%t}(-1 to -!k)
				equation eq!kb_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-6)) fs_!j1_{%y}_{%t}(-6 to !k6) fs_!j2_{%y}_{%t}(-6 to !k6)
				equation eq!kc_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-12)) fs_!j1_{%y}_{%t}(-12 to !k12) fs_!j2_{%y}_{%t}(-12 to !k12)
			else
			if !j = 3 then
				equation eq!ka_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-1)) fs_!j1_{%y}_{%t}(-1 to -!k) fs_!j2_{%y}_{%t}(-1 to -!k) fs_!j3_{%y}_{%t}(-1 to -!k)
				equation eq!kb_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-6)) fs_!j1_{%y}_{%t}(-6 to !k6) fs_!j2_{%y}_{%t}(-6 to !k6) fs_!j3_{%y}_{%t}(-6 to !k6)
				equation eq!kc_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-12)) fs_!j1_{%y}_{%t}(-12 to !k12) fs_!j2_{%y}_{%t}(-12 to !k12) fs_!j3_{%y}_{%t}(-12 to !k12)
			else
			if !j = 4 then
				equation eq!ka_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-1)) fs_!j1_{%y}_{%t}(-1 to -!k) fs_!j2_{%y}_{%t}(-1 to -!k) fs_!j3_{%y}_{%t}(-1 to -!k) fs_!j4_{%y}_{%t}(-1 to -!k)
				equation eq!kb_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-6)) fs_!j1_{%y}_{%t}(-6 to !k6) fs_!j2_{%y}_{%t}(-6 to !k6) fs_!j3_{%y}_{%t}(-6 to !k6) fs_!j4_{%y}_{%t}(-6 to !k6)
				equation eq!kc_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-12)) fs_!j1_{%y}_{%t}(-12 to !k12) fs_!j2_{%y}_{%t}(-12 to !k12) fs_!j3_{%y}_{%t}(-12 to !k12) fs_!j4_{%y}_{%t}(-12 to !k12)
			else
			if !j = 5 then
				equation eq!ka_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-1)) fs_!j1_{%y}_{%t}(-1 to -!k) fs_!j2_{%y}_{%t}(-1 to -!k) fs_!j3_{%y}_{%t}(-1 to -!k) fs_!j4_{%y}_{%t}(-1 to -!k) fs_!j5_{%y}_{%t}(-1 to -!k)
				equation eq!kb_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-6)) fs_!j1_{%y}_{%t}(-6 to !k6) fs_!j2_{%y}_{%t}(-6 to !k6) fs_!j3_{%y}_{%t}(-6 to !k6) fs_!j4_{%y}_{%t}(-6 to !k6) fs_!j5_{%y}_{%t}(-6 to !k6)
				equation eq!kc_!j_{%y}_{%t}.ls(cov=hac) @pcy(us) @pcy(us(-12)) fs_!j1_{%y}_{%t}(-12 to !k12) fs_!j2_{%y}_{%t}(-12 to !k12) fs_!j3_{%y}_{%t}(-12 to !k12) fs_!j4_{%y}_{%t}(-12 to !k12) fs_!j5_{%y}_{%t}(-12 to !k12)
			endif 
			endif
			endif
			endif
			next	
		next
	next
next
'Estimate equations over various specifications and lags


'%%% FORECAST USING ESTIMATED EQUATIONS

delete fc*
' Get rid of old forecasts

for %y reg all
	for %t cov cor
		for !j = 2 to 5
			for !k = 1 to !max_lag
				for !h = 0 to 11
					smpl 2012m06+!h 2012m06+!h

					eq!ka_!j_{%y}_{%t}.forecast(g,e) fc!ka_!j_{%y}_{%t}_t
					eq!kb_!j_{%y}_{%t}.forecast(g,e) fc!kb_!j_{%y}_{%t}_t
					eq!kc_!j_{%y}_{%t}.forecast(g,e) fc!kc_!j_{%y}_{%t}_t

					genr fc!ka_!j_{%y}_{%t} = fc!ka_!j_{%y}_{%t}_t
					genr fc!kb_!j_{%y}_{%t} = fc!kb_!j_{%y}_{%t}_t
					genr fc!kc_!j_{%y}_{%t} = fc!kc_!j_{%y}_{%t}_t

					delete *_t
				next
			next
		next
	next
next
' Generate monthly rolling forecasts over all specifications at 1,6, and 12 month horizons.

!num_rows = !max_lag *16+1
table(!num_rows,10) results_table
' This has 2 different aggregations times maf from 2 different lags times 4 different factor inclusions times 3 different lag settings for forecast equations.

results_table(1,1) = "Specification"
results_table(1,2) = "RMSE_1m"
results_table(1,3) = "MAE_1m"
results_table(1,4) = "MAPE_1m"
results_table(1,5) = "RMSE_6m"
results_table(1,6) = "MAE_6m"
results_table(1,7) = "MAPE_6m"
results_table(1,8) = "RMSE_12m"
results_table(1,9) = "MAE_12m"
results_table(1,10) = "MAPE_12m"
' Table headings for columns

smpl 2012m06 2013m05
' Evaluate forecasts for one year forward

!counter = 2
for %y reg all
	for %t cov cor
		for !k = 1 to !max_lag
			for !j = 2 to 5
				%j = @str(!j)
				%k = @str(!k)

				results_table(!counter,1) = "fc" + %k + "_"+ %j + "_" + %y + "_" + %t
				results_table(!counter,2) = @rmse(fc!ka_!j_{%y}_{%t},us)
				results_table(!counter,3) = @mae(fc!ka_!j_{%y}_{%t},us)
				results_table(!counter,4) = @mape(fc!ka_!j_{%y}_{%t},us)
				results_table(!counter,5) = @rmse(fc!kb_!j_{%y}_{%t},us)
				results_table(!counter,6) = @mae(fc!kb_!j_{%y}_{%t},us)
				results_table(!counter,7) = @mape(fc!kb_!j_{%y}_{%t},us)
				results_table(!counter,8) = @rmse(fc!kc_!j_{%y}_{%t},us)
				results_table(!counter,9) = @mae(fc!kc_!j_{%y}_{%t},us)
				results_table(!counter,10) = @mape(fc!kc_!j_{%y}_{%t},us)

				!counter = !counter+1
			next
		next
	next
next

results_table.save(t=csv, n="NA") dfm_g2_fc_results
' Save results and export to excel

smpl 2007m01 @last