# Yet Another Blog in Statistical Computing

I can calculate the motion of heavenly bodies but not the madness of people. -Isaac Newton

## A SAS Macro for Breusch-Pagan Test

In SAS, Breusch-Pagan test for Heteroscedasticity in a linear regression can be conducted with MODEL procedure in SAS/ETS, as shown in the code snippet below.

```data one;
do i = 1 to 100;
x1 = uniform(1);
x2 = uniform(2);
r  = normal(1) * 0.1;
if x2 > 0.5 then r = r * 2;
y = 10 + 8 * x1 + 6 * x2 + r;
output;
end;
run;

proc model data = one;
parms b0 b1 b2;
y = b0 + b1 * x1 + b2 * x2;
fit y / breusch = (1 x1 x2);
run;
/*
Heteroscedasticity Test
Equation        Test               Statistic     DF    Pr > ChiSq    Variables

y               Breusch-Pagan          10.44      2        0.0054    1, x1, x2
*/
```

However, in a forecasting model that I am recently working on, I find that it is not convenient to use “proc model” every time when I want to do Breusch-Pagan test and rather prefer a more generic solution not tied to a specific SAS module or procedure and that would only need to take a minimum set of inputs instead of specifying out a full model. As a result, I draft a simple sas macro to do Breusch-Pagan test, which gives the identical result as the one from MODEL procedure. Hopeful, others might find this macro useful as well.

```proc reg data = one;
model y = x1 x2;
output out = two r = resid;
run;

%macro hetero_bp(r = , x = , data = );
***********************************************************;
* THE SAS MACRO IS TO CALCULATE BREUSCH-PAGEN TEST FOR    *
* HETEROSKEDASTICITY                                      *;
* ======================================================= *;
* PAMAMETERS:                                             *;
*  DATA: INPUT SAS DATA TABLE                             *;
*  R   : RESIDUAL VALUES FROM A LINEAR REGRESSION         *;
*  X   : A LIST OF NUMERIC VARIABLES TO MODEL ERROR       *;
*        VARIANCE IN BREUSCH-PAGEN TEST                   *;
* ======================================================= *;
* CONTACT:                                                *;
*  WENSUI.LIU@53.COM                                      *;
***********************************************************;

data _data_(keep = r2 &x);
set &data;
where r ~= .;
r2 = &r ** 2;
run;

ods output nobs = _nobs_;
ods output anova = _anova_;
ods output fitstatistics = _fits_;
ods listing close;
proc reg data = _last_;
model r2 = &x;
run;
ods listing;

proc sql noprint;
select distinct NObsUsed into :n from _nobs_;
select df into :df from _anova_ where upcase(compress(source, ' ')) = 'MODEL';
select nvalue2 into :r2 from _fits_ where upcase(compress(label2, ' ')) = 'R-SQUARE';
run;
%put &r2;

data _result_;
chi_square = &r2 * &n;
df         = &df;
p_value    = 1 - probchi(chi_square, df);
run;

proc report data = _last_ spacing = 1 headline nowindows split = "*";
column(" * BREUSCH-PAGEN TEST FOR HETEROSKEDASTICITY * "
chi_square df p_value);
define chi_square  / "CHI-SQUARE"  width = 15;
define df          / "DF"          width = 5;
define p_value     / "P-VALUE"     width = 15 format = 12.8;
run;

proc datasets library = work;
delete _: / memtype = data;
run;

%mend hetero_bp;

%hetero_bp(r = resid, x = x1 x2, data = two);
/*
BREUSCH-PAGEN TEST FOR HETEROSKEDASTICITY

CHI-SQUARE    DF         P-VALUE
-----------------------------------------
10.4389     2      0.00541030
*/
```