/* Simulate various types of time series data */
data sim;
call streaminit(12345);
do t=1 to 5000;
	e = rand("Normal",0,1);
	u = rand("Normal",0,1);
	x = rand("Normal",2,2);

	/* Serial autocorrelation, AR(1) */
	if t=1 then do;
		e = u;					* Autocorrelation structure ;
		y = 1 + e;				* Generate data ;
	end;
	else do;
		e = 0.9*e_l + u;		* Autocorrelation structure ;
		y = 1 + 1*x + e;		* Generate data ; 
	end;
	output;

	/* Generate lagged values of errors */
	e_l = e;

end;
label y = "Autocorrelated data"
	  t = "Time";
run;

/* Draw a random sample */
proc surveyselect data=sim
  method=srs n=500 seed=123
  out=sim_sample noprint;
run;

/* Testing for serial autocorrelation -- Durbin-Watson test*/
proc reg data=sim_sample;
model y = x / dwprob;
run;

/* Correcting for autocorrelation */
/* (1) FGLS */

	/* Initial regression */
	proc reg data=sim_sample plots=(none);
	model y=x;
	output out=reg_e r=e_hat;
	ods output ParameterEstimates=sim_sac_se(keep=estimate stderr);
	run;quit;

	/* Determine DW rho_hat */
	data reg_e;
	set reg_e;
	e_hatl = lag(e_hat);
	eel = e_hat * e_hatl;
	ee = e_hat * e_hat;
	run;
	proc means data=reg_e noprint;
	var ee eel;
	output out=rho_hat(drop=_TYPE_ _FREQ_) sum=ee eel;
	run;
	data sim_correct;
	set sim_sample;
	if _n_ = 1 then set rho_hat;
	rho_hat = ((500**2)*(eel/ee)/2+(2**2))/((500**2)-(2**2));
	yl = rho_hat*lag(y);
	xl = rho_hat*lag(x);
	y_star = y - yl;
	x_star = x - xl;
	int_star = 1 - rho_hat;
	run;

	/* Estimate corrected OLS regression */
	proc reg data=sim_correct plots=(none);
	model y_star = int_star x_star / noint dwprob;
	ods output ParameterEstimates=sim_fgls_se(keep=estimate stderr rename=(estimate=estimate_fgls stderr=stderr_fgls))
	run;
	
/* (2) Newey-West robust standard errors */
	proc model data=sim_sample;
    endo y;
    exog x;
    instruments _exog_;
    parms b0 b1;
    y = b0 + b1*x;
    fit y / gmm kernel=(bart,1,0) vardef=n dw; 	* Second parameter in kernel specifies lag length ; 
	ods output ParameterEstimates=sim_nw_se(keep=estimate stderr rename=(estimate=estimate_nw stderr=stderr_nw));
    run;quit; 

/* View and compare standard errors from each regression */
data sim_se;
merge sim_sac_se sim_fgls_se sim_nw_se;
label estimate = "Parameters (OLS)"
	  stderr = "SE (OLS)"
	  estimate_fgls = "Parameters (FGLS)"
	  stderr_fgls = "SE (FGLS)"
	  estimate_nw = "Parameters (Newey-West)"
	  stderr_nw = "SE (Newey-West)";
run;
proc print data=sim_se noobs label;
run;