/* Summary statistics */
proc means data=rma_canola mean median std min max;
var cov_lvl policies_sold net_acres total_prem liability indem_amt loss_ratio;
run;

/* Summation across insurance policies */
proc sort data=rma_canola;
by fips cov_lvl;
run;
proc means data=rma_canola noprint;
by fips cov_lvl;
output out=rma_canola(drop=_TYPE_ _FREQ_) sum=;
run;

/* Eliminate observations where policies>1 but net acres insured = 0;
	*/
data rma_canola;
set rma_canola;
if policies_sold>0 and net_acres=0 then delete;
run;

/* OLS regression */
proc reg data=rma_canola;
model policies_sold = net_acres total_prem liability loss_ratio;
run;quit;

/* OLS regression, residual analysis */
ods graphics on;
proc reg data=rma_canola plots=all;
model policies_sold = net_acres total_prem liability loss_ratio;
run;quit;
ods graphics off;

/* Test for heteroskedasticity; Null is homoskedasticity */
proc reg data=rma_canola;
model policies_sold = net_acres total_prem liability 
					  loss_ratio / spec ;
run;quit;

/* Perform FGLS */
/* Generate residuals */
proc reg data=rma_canola noprint;
model policies_sold = net_acres total_prem liability 
					  loss_ratio ;
output out=resids r=e_hat;
run;quit;

/* Log squared residuals */
data resids;
set resids;
e_hat = log(e_hat**2);
run;

/* Regress logged residuals on X;
	generate predicted values */
proc reg data=resids noprint;
model e_hat = net_acres total_prem liability 
					  loss_ratio ;
output out=predict p=e_pred;
run;quit;

/* Create weights */
data rma_canola_weighted;
set predict;
weights = 1/sqrt(e_pred);
run;

/* Run a weighted regression */
proc reg data=rma_canola_weighted;
weight weights;
model policies_sold = net_acres total_prem liability 
					  loss_ratio ;
run;quit;

/* OLS with White's heteroskedasiticity robust errors */
proc reg data=rma_canola;
model policies_sold = net_acres total_prem liability 
					  loss_ratio/white hccmethod=2 ;
run;quit;