workfile DTR 7 1/1/1960  7/31/2018
	
	smpl 1/1/1960 7/31/2018

	
	'Create trend and monthly dummies
	
	genr time = @trend + 1

	genr d1=@month=1
	genr d2=@month=2
	genr d3=@month=3
	genr d4=@month=4
	genr d5=@month=5
	genr d6=@month=6
	genr d7=@month=7
	genr d8=@month=8
	genr d9=@month=9
	genr d10=@month=10
	genr d11=@month=11
	genr d12=@month=12
	
	
	for %station atl bos bwi cvg dfw dsm dtw  las lga  msp ord pdx  phl slc tus

	
	'Read data and create AVG and DTR
	
	read {%station}.txt max{%station} min{%station}
	
	genr AVG_{%station}=(max{%station}+min{%station})/2
	genr DTR_{%station}=max{%station}-min{%station}

	
	'Kernel density estimates
	
	smpl 1960 2017
	
	freeze(a0_{%station}) AVG_{%station}.distplot kernel
	a0_{%station}.setelem(1)  lwidth(3)
	
	freeze(a00_{%station}) DTR_{%station}.distplot kernel
	a00_{%station}.setelem(1)  lwidth(3)
	
	graph a000_{%station}.merge a0_{%station} a00_{%station}
	a000_{%station}.align(2,1.5,0) 
	a000_{%station}.axis(bottom) range(-10,100)
	a000_{%station}.axis(left) range(0,.065)

	
	'***** Sequential regressions *****
	
	
	'Trend 
	
	smpl 1960 2017
	
	equation a155_{%station}
	
	a155_{%station}.ls(cov=hac) AVG_{%station} c time  
	freeze(taba155_{%station}) a155_{%station}.stats
	
	genr AVGtrend_{%station} = AVG_{%station} - resid
	
	group groupxxx_{%station} AVG_{%station} AVGtrend_{%station}
	freeze(a1_{%station}) groupxxx_{%station}.line
	a1_{%station}.setelem(1) lcolor(gray)
	a1_{%station}.setelem(2)  lwidth(3)  lcolor(blue)
	a1_{%station}.legend -display
	
	genr AVGdet_{%station} = resid
	
	equation a156_{%station}
	
	a156_{%station}.ls(cov=hac)  DTR_{%station} c time 
	freeze(taba156_{%station}) a156_{%station}.stats
	
	genr DTRtrend_{%station} = DTR_{%station} - resid
	
	group groupyyy_{%station} DTR_{%station} DTRtrend_{%station}
	freeze(a11_{%station}) groupyyy_{%station}.line
	a11_{%station}.setelem(1) lcolor(gray)
	a11_{%station}.setelem(2)  lwidth(3) lcolor(blue)
	a11_{%station}.legend -display
	
	genr DTRdet_{%station} = resid
	
	graph a111_{%station}.merge a1_{%station} a11_{%station}
	a111_{%station}.align(2,1.5,0) 
	
	
	'Fixed seasonality
	
	smpl 1960 2017
	
	equation a444a_{%station}
	
	a444a_{%station}.ls(cov=hac) AVGdet_{%station} d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12
	freeze(taba444a_{%station}) a444a_{%station}.stats
	
	genr AVGdetseas_{%station} = AVGdet_{%station} - resid
	
	group groupy2a_{%station} avgdet_{%station} avgdetseas_{%station}
	freeze(a4444a_{%station}) groupy2a_{%station}.line
	a4444a_{%station}.setelem(1) lcolor(gray)
	a4444a_{%station}.setelem(2)  lwidth(3) lcolor(blue)
	a4444a_{%station}.options size(20,3)
	a4444a_{%station}.legend -display
	
	genr AVGdetdes_{%station}=resid
	
	equation a444_{%station}
	a444_{%station}.ls(cov=hac) DTRdet_{%station} d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12
	
	freeze(taba444_{%station}) a444_{%station}.stats
	
	genr DTRdetseas_{%station} = DTRdet_{%station} - resid
	
	group groupy2_{%station} DTRdet_{%station} DTRdetseas_{%station}
	freeze(a4444_{%station}) groupy2_{%station}.line
	a4444_{%station}.setelem(1) lcolor(gray)
	a4444_{%station}.setelem(2)  lwidth(3) lcolor(blue)
	a4444_{%station}.options size(20,3)
	a4444_{%station}.legend -display
	
	genr DTRdetdes_{%station}=resid
	
	graph a4444b_{%station}.merge a4444a_{%station} a4444_{%station}

	
	'12-month fixed seasonal patterns
	
	smpl 2017 2017
	
	freeze(a6666db_{%station}) AVGdetseas_{%station}.line
	a6666db_{%station}.setelem(1)  lwidth(3)
	
	freeze(a6666d_{%station}) DTRdetseas_{%station}.line
	a6666d_{%station}.setelem(1)  lwidth(3)
	
	graph a6666ddd_{%station}.merge a6666db_{%station} a6666d_{%station}
	a6666ddd_{%station}.align(2,1.5,0) 
		
	
	'Time-varying seasonality: full-sample fits
	
	smpl 1960 2017
	
	equation a666AVG_{%station}
	
	a666AVG_{%station}.ls(cov=hac) AVGdet_{%station} d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d1*time d2*time  d3*time  d4*time  d5*time  d6*time  d7*time  d8*time  d9*time  d10*time  d11*time  d12*time 
	
	freeze(taba666AVG_{%station}) a666AVG_{%station}.stats
	
	genr AVGdetseas2_{%station} = AVGdet_{%station} - resid
	genr AVGdetdes2_{%station} = resid
		
	equation a666_{%station}
	
	a666_{%station}.ls(cov=hac) DTRdet_{%station} d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d1*time d2*time  d3*time  d4*time  d5*time  d6*time  d7*time  d8*time  d9*time  d10*time  d11*time  d12*time 
	
	freeze(taba666_{%station}) a666_{%station}.stats
	
	genr DTRdetseas2_{%station} = DTRdet_{%station} - resid
	genr DTRdetdes2_{%station} = resid
		
	
	'Time-Varying Seasonality:  Fits at beginning and end
	
	smpl 1/1/2017 12/31/2017
	
	genr earlyshiftAVG_{%station} = AVGdetseas2_{%station}(-20820)
	
	group patternsAVG_{%station} earlyshiftAVG_{%station} AVGdetseas2_{%station}
	freeze(help2AVG_{%station}) patternsAVG_{%station}.line
	help2AVG_{%station}.setelem(1)  lwidth(3)  lcolor(blue)
	help2AVG_{%station}.setelem(2)  lwidth(3)  lcolor(red)
	help2AVG_{%station}.legend -display
	
	
	genr earlyshiftDTR_{%station} = DTRdetseas2_{%station}(-20820)
	
	group patternsDTR_{%station} earlyshiftDTR_{%station} DTRdetseas2_{%station}
	freeze(help2DTR_{%station}) patternsDTR_{%station}.line
	help2DTR_{%station}.setelem(1)  lwidth(3)  lcolor(blue)
	help2DTR_{%station}.setelem(2)  lwidth(3)  lcolor(red)
	help2DTR_{%station}.legend -display
	
	graph help2222_{%station}.merge help2AVG_{%station} help2DTR_{%station}
	help2222_{%station}.align(2,1.5,0) 
	
		
	'***** JOINT regressions *****
	
	smpl 1960 2017
	
	'Conditional Mean Regression, AVG
	
	equation avg_np1_{%station}
	avg_np1_{%station}.ls(cov=hac) avg_{%station}  c time  d1 d2 d3 d4 d5 d6   d8 d9 d10 d11 d12 d1*time  d2*time  d3*time  d4*time  d5*time  d6*time   d8*time  d9*time  d10*time  d11*time  d12*time avg_{%station}(-1) 
	
	genr e_avg_np1_{%station}=resid
	
	'p(nt)
	freeze(avg_np1_nt_{%station}) avg_np1_{%station}.wald c(2)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(ns)
	freeze(avg_np1_ns_{%station}) avg_np1_{%station}.wald c(3)=c(4)=c(5)=c(6)=c(7)=c(8)=c(9)=c(10)=c(11)=c(12)=c(13)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(nts)
	freeze(avg_np1_nts_{%station}) avg_np1_{%station}.wald c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0 
	
	freeze(taba777cc_{%station}) avg_np1_{%station}.stats
	
	
	'Conitional Mean Regression, DTR
	
	equation dtr_np1_{%station}
	dtr_np1_{%station}.ls(cov=hac) dtr_{%station} c time d1 d2 d3 d4 d5 d6  d8 d9 d10 d11 d12 d1*time  d2*time  d3*time  d4*time  d5*time  d6*time   d8*time  d9*time  d10*time  d11*time  d12*time  dtr_{%station}(-1)
	
	genr e_dtr_np1_{%station}=resid
	
	'p(nt)
	freeze(dtr_np1_nt_{%station}) dtr_np1_{%station}.wald c(2)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(ns)
	freeze(dtr_np1_ns_{%station}) dtr_np1_{%station}.wald c(3)=c(4)=c(5)=c(6)=c(7)=c(8)=c(9)=c(10)=c(11)=c(12)=c(13)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(nts)
	freeze(dtr_np1_nts_{%station}) dtr_np1_{%station}.wald c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0 
	
	freeze(taba777dd_{%station}) dtr_np1_{%station}.stats
		
	
	'Conditional Variance Regression, AVG
	
	genr e2_avg_{%station} = (e_avg_np1_{%station})^2
	
	equation avg_np2_{%station}
	avg_np2_{%station}.ls(cov=hac) e2_avg_{%station}  c time d1 d2 d3 d4 d5 d6   d8 d9 d10 d11 d12 d1*time  d2*time  d3*time  d4*time  d5*time  d6*time   d8*time  d9*time  d10*time  d11*time  d12*time e2_avg_{%station}(-1) 
	
	'p(nt)
	freeze(avg_np2_nt_{%station}) avg_np2_{%station}.wald c(2)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(ns)
	freeze(avg_np2_ns_{%station}) avg_np2_{%station}.wald c(3)=c(4)=c(5)=c(6)=c(7)=c(8)=c(9)=c(10)=c(11)=c(12)=c(13)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(nts)
	freeze(avg_np2_nts_{%station}) avg_np2_{%station}.wald c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0 
	
	freeze(taba777cccc_{%station}) avg_np2_{%station}.stats
	
	
	'Conditional Variance Regression, DTR
	
	genr e2_dtr_{%station} = (e_dtr_np1_{%station})^2
	
	equation dtr_np2_{%station}
	dtr_np2_{%station}.ls(cov=hac) e2_dtr_{%station}  c time d1 d2 d3 d4 d5 d6   d8 d9 d10 d11 d12 d1*time  d2*time  d3*time  d4*time  d5*time  d6*time   d8*time  d9*time  d10*time  d11*time  d12*time e2_dtr_{%station}(-1) 
	
	'p(nt)
	freeze(dtr_np2_nt_{%station}) dtr_np2_{%station}.wald c(2)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(ns)
	freeze(dtr_np2_ns_{%station}) dtr_np2_{%station}.wald c(3)=c(4)=c(5)=c(6)=c(7)=c(8)=c(9)=c(10)=c(11)=c(12)=c(13)=c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0
	
	'p(nts)
	freeze(dtr_np2_nts_{%station}) dtr_np2_{%station}.wald c(14)=c(15)=c(16)=c(17)=c(18)=c(19)=c(20)=c(21)=c(22)=c(23)=c(24)=0 
	
	freeze(taba777dddd_{%station}) dtr_np2_{%station}.stats
		
	
	'Standardized residual skewness and kurtosis
	
	avg_np2_{%station}.fit sigma2_avg_{%station}
	genr sigma_avg_{%station} = sigma2_avg_{%station}^.5
	genr  standardized_avg_{%station} =  e_avg_np1_{%station} / sigma_avg_{%station}
	
	dtr_np2_{%station}.fit sigma2_dtr_{%station}
	genr sigma_dtr_{%station} = sigma2_dtr_{%station}^.5
	genr  standardized_dtr_{%station} =  e_dtr_np1_{%station} / sigma_dtr_{%station}
	
	next
	
	stop