Repository: JuliaStats/HypothesisTests.jl
Pull request #206 (2020-05-26 14:24:38)

add White's test (heteroskedasticity)
This PR adds White's (and Breusch-Pagan's) test for heteroskedasticity in a regression. It is one of the most commonly applied diagnostic tests (along with various tests for autocorrelation which already are in the package). The test gives the same result as `bptest` in R.

Closes issue https://github.com/JuliaStats/HypothesisTests.jl/issues/199

azev77: Thanks this looks great! It might be nice to explicitly define `BreuschPaganTest` in case someone looking for it doesn't realize it's a special case of White Test. 
Something that reuses code like: 
```julia
BreuschPaganTest(X, e) = WhiteTest(X, e, :Linear)
```
PaulSoderlind: @azev77 Thanks for the feedback. One way of doing that is to 

(1) add the function as you suggest, export it
(2) output a string from `WhiteTest()` that describes the test (eg. `TestDescription = "Breusch-Pagan's test for heteroskedasticity (linear)"` depending on the `TestType`
(3) which is used for `testname(t::WhiteTest) = t.TestDescription`

Sounds OK?
azev77: Sounds good
PaulSoderlind: added Breusch-Pagan test and more precise `testname(t::WhiteTest)`
azev77: Is ‘ population_param_of_interest’ standard? 
I like ‘estimand’. 
I’m not recommending a modification, just discussing for the future...
PaulSoderlind: @azev77 Thanks. I just tried to follow what appears to be the default in this package.

Another thing: it would be nice to make this package and `CovarianceMatrices.jl` communicate better. The easiest way is to make something like https://github.com/gragusa/CovarianceMatrices.jl/issues/49 happen. Please drop me a line in case you are interested in contributing.
azev77: @PaulSoderlind 
check out a [recent paper](https://www.sciencedirect.com/science/article/abs/pii/S030440761630197X) (Romano & Wolf, JE, 2017)
"Therefore, it is tempting to decide based on the data which route to take: OLS or WLS. Specifically,we suggest applying a test for conditional heteroskedasticity.  Several such tests exists, the mostpopular ones being the tests of **Breusch and Pagan** (1979) and **White** (1980); also see **Koenker**(1981) and **Koenker and Bassett** (1982).  If the null hypothesis of conditional homoskedasticity itnot rejected by such a test, use the OLS estimator; otherwise, use the WLS estimator.  We call theresulting estimator theadaptive least squares(ALS) estimator.  The motivation is as follows.  Underconditional homoskedasticity, the ALS estimator will be equal to the WLS estimator with a smallprobability only (roughly equal to the nominal size of the test).  Therefore, in this case, ALS isexpected to be more efficient than WLS in finite samples, though still less efficient than OLS."
azev77: @PaulSoderlind 
in your code you regress `e2` on `z=[X, X.^2]`
I've [seen some](https://jsf.pm-research.com/content/22/1/35.abstract) regress `e2` on `z=[y_hat, y_hat^2`
PaulSoderlind: @azev77 Yes, there are lots of versions of these tests. For instance, what I called the BP test is really the BP/Koenker (1981) test.

The code above can be used to perform some of the other versions applied in other packages: just specify `X` as desired to have 
1.  just some of the original regressors 
2. the predicted values (`y_hat`) instead of the original regressors. 

Maybe we should add this information to the documentation of the function?

As you point out, these are restricted versions, sometimes applied save dfs/increase power. I believe Stata has some interesting routines for this, and even through on a Bonferroni approach for multiple testing. I have Julia code for Bonferroni/Holm so that could be done here - at some point in time.
PaulSoderlind: @wildart I added a comment along the lines of the comment by @azev77.

Good to go?
wildart: You forgot to add the test file to `test/runtests.jl`
PaulSoderlind: @wildart I added white.jl to `test/runtests.jl`. I also added `BreuschPaganTest` to doc and test.

The test passes. I did notice, however, that `Pkg.test(HypothesisTests)` fails because of `show.jl`. This appears to be unrelated to my PR.
wildart: Disregard nightly builds fails, especially if it doesn't related to your PR.
PaulSoderlind: It seems the feedback so far is positive. Anything else I can do to make this move forward?
wildart: @ararslan @mschauer Can we merge this PR?
PaulSoderlind: @nalimilan 

Here is the problem with keyword arguments that I referred to above (my setting is Win 10, Julia 1.4.2/1.5.0-beta1.0)

```
struct foo
    A::Int
    B::Symbol
end

#function foo(x,B=:bar)          #works
function foo(x;B=:bar)           #gives InexactError: Int64(1.2)
  x = round(Int,x)
  return foo(x,B)
end

z = foo(1.2,:xx)
```
nalimilan: Sorry, I don't understand the problem. You need to write `z = foo(1.2, B=:xx)` when using keyword arguments. `z = foo(1.2, :xx)` calls the default constructor, which you can override if needed using an inner constructor (see [the manual](https://docs.julialang.org/en/latest/manual/constructors/#man-inner-constructor-methods-1)). But I don't think that's the case here?
PaulSoderlind: @nalimilan I have made the changes you requested and make an update of the PR. And, yes, you are right about the keyword argument. 
PaulSoderlind: @nalimilan and @azev77 

I have some **2nd thoughts** about one of latest design changes: `X` should really by an `AbstractMatrix` (not an `AbstractArray`) and we should check whether it actually includes a (non-zero) constant column. Motivation: the test really needs `X` to have a constant and at least one more column. If not, the results are not valid (so a silent failure). 

To test the `X` matrix: (comments welcome):
```
constant_cols = mapslices(col->all(diff(col).==0),X,dims=1)   #loops over columns, test if constant
if !any(X[1,vec(constant_cols)] .!= 0)          #test if any constant column is non-zero
  throw(ArgumentError("one of the colums of X must be a non-zero constant"))
end
```

This clearly allocates memory, but just one column at a time.
PaulSoderlind: I have now implemented (almost all of) the suggestions by @nalimilan limilan. In addition, I have added checks on `X` having at least two columns and including a non-zero constant column (or else the test can silently fail). Tests for these checks have been added.

@wildart good to go? 
PaulSoderlind: @nalimilan Thanks for the suggestions. I believe I implemented all of them. Yes, your `intercept_col` is a better idea than my old code. I changed slightly to
`intercept_col = any(col -> first(col) != 0 && all(==(first(col)), col),eachcol(X))`

However, I now notice travis test failures on Julia 1.0 (passes the 1.3 test). Because of `eachcol`?

nalimilan: > ``
> 
> However, I now notice travis test failures on Julia 1.0 (passes the 1.3 test). Because of `eachcol`?

Ah, yes, you'll need to add a dependency on Compat 2.2.0.

Patch:
diff --git a/docs/src/time_series.md b/docs/src/time_series.md
index e3f8e96..ab3865a 100644
--- a/docs/src/time_series.md
+++ b/docs/src/time_series.md
@@ -27,4 +27,10 @@ ADFTest
 ## Diebold-Mariano test
 ```@docs
 DieboldMarianoTest
-```
\ No newline at end of file
+```
+
+## White test
+```@docs
+WhiteTest
+BreuschPaganTest
+```
diff --git a/src/white.jl b/src/white.jl
new file mode 100644
index 0000000..93ba8a9
--- /dev/null
+++ b/src/white.jl
@@ -0,0 +1,120 @@
+# white.jl
+# White and Breusch-Pagan tests for heteroskedasticity
+#
+# Copyright (C) 2020   Paul Söderlind
+#
+# Permission is hereby granted, free of charge, to any person obtaining
+# a copy of this software and associated documentation files (the
+# "Software"), to deal in the Software without restriction, including
+# without limitation the rights to use, copy, modify, merge, publish,
+# distribute, sublicense, and/or sell copies of the Software, and to
+# permit persons to whom the Software is furnished to do so, subject to
+# the following conditions:
+#
+# The above copyright notice and this permission notice shall be
+# included in all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+export WhiteTest, BreuschPaganTest
+
+struct WhiteTest <: HypothesisTest
+    dof::Int                 #degrees of freedom in test
+    lm::Float64              #Lagrange Multiplier test statistic, distributed as Chisq(dof)
+    type::Symbol
+end
+
+"""
+    WhiteTest(X, e; type = :White)
+
+Compute White's (or Breusch-Pagan's) test for heteroskedasticity.
+
+`X` is a matrix of regressors and `e` is the vector of residuals from the original model.
+The keyword argument `type` is either `:linear` for the Breusch-Pagan/Koenker test,
+`:linear_and_squares` for White's test with linear and squared terms only (no cross-products), or
+`:White` (the default) for the full White's test (linear, squared and cross-product terms).
+`X` should include a constant and at least one more regressor, with observations in rows and
+regressors in columns. In some applications, `X` is a subset of the regressors in
+the original model, or just the fitted values. This saves degrees of freedom and may
+give a more powerful test. The `lm` (Lagrange multiplier) test statistic is T*R2
+where R2 is from the regression of `e^2` on the terms mentioned above.
+Under the null hypothesis it is distributed as `Chisq(dof)` where `dof` is the number
+of independent terms (not counting the constant), so the null is rejected when the
+test statistic is large enough.
+
+Implements: [`pvalue`](@ref)
+
+# References
+  * H. White, (1980): A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity, Econometrica, 48, 817-838.
+  * T.S. Breusch & A.R. Pagan (1979), A simple test for heteroscedasticity and random coefficient variation, Econometrica, 47, 1287-1294
+  * R. Koenker (1981), A note on studentizing a test for heteroscedasticity, Journal of Econometrics, 17, 107-112
+
+# External links
+  * [White's test on Wikipedia](https://en.wikipedia.org/wiki/White_test)
+"""
+function WhiteTest(X::AbstractMatrix{<:Real}, e::AbstractVector{<:Real}; type = :White)
+    (n,K) = (size(X,1),size(X,2))    #nobs,nvars
+    n == length(e) || throw(DimensionMismatch("inputs must have the same length"))
+
+    K >= 2 || throw(ArgumentError("X must have >= 2 columns"))
+
+    intercept_col = false
+    for i = 1:K
+        col           = view(X,:,i)
+        intercept_col = first(col) != 0 && all(==(first(col)), col)
+        intercept_col && break
+    end
+    intercept_col || throw(ArgumentError("one of the colums of X must be a non-zero constant"))
+
+    if type == :linear                    #Breush-Pagan/Koenker
+        z = X
+    elseif type == :linear_and_squares    #White with linear and squares
+        z = [X X.^2]
+    else                                      #White with linear, squares and cross-products
+        z = fill(NaN,n,round(Int,K*(K+1)/2))
+        vv = 1
+        @views for i = 1:K, j = i:K
+            z[:,vv] .= X[:,i] .* X[:,j]       #eg. x1*x1, x1*x2, x2*x2
+            vv      += 1
+        end
+    end
+
+    dof = rank(z) - 1                        #number of independent regressors in z
+    e2  = e.^2
+    b   = z\e2
+    res = e2 - z*b
+    R2  = 1 - var(res)/var(e2)
+    lm  = n*R2
+    return WhiteTest(dof, lm, type)
+end
+
+"""
+    BreuschPaganTest(X, e)
+
+Compute Breusch-Pagan's test for heteroskedasticity.
+
+`X` is a matrix of regressors from the original model and `e` the vector of residuals.
+This is equivalent to `WhiteTest(X, e, type = :linear)`.
+See [`WhiteTest`](@ref) for further details.
+"""
+BreuschPaganTest(X, e) = WhiteTest(X, e, type = :linear)
+
+testname(t::WhiteTest) = "White's (or Breusch-Pagan's) test for heteroskedasticity"
+
+population_param_of_interest(t::WhiteTest) = ("T*R2", 0, t.lm)
+default_tail(test::WhiteTest)              = :right
+
+function show_params(io::IO, t::WhiteTest, ident = "")
+    println(io, ident, "T*R^2 statistic:        ", t.lm)
+    println(io, ident, "degrees of freedom:     ", t.dof)
+    println(io, ident, "type:                   ", t.type)
+end
+
+StatsBase.dof(t::WhiteTest) = t.dof
+pvalue(t::WhiteTest) = pvalue(Chisq(t.dof), t.lm, tail=:right)


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index 4a5cf22..c1b5660 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -191,5 +191,6 @@ include("wald_wolfowitz.jl")
 include("f.jl")
 include("correlation.jl")
 include("diebold_mariano.jl")
+include("white.jl")
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index a2f91a4..cf06c82 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -37,3 +37,4 @@ include("wilcoxon.jl")
 include("z.jl")
 include("f.jl")
 include("diebold_mariano.jl")
+include("white.jl")
diff --git a/test/white.jl b/test/white.jl
new file mode 100644
index 0000000..df11865
--- /dev/null
+++ b/test/white.jl
@@ -0,0 +1,252 @@
+using HypothesisTests, Test
+
+@testset "White tests" begin
+
+X = [ 4.18      0.77      1.00;
+     -3.41      0.73      1.00;
+      5.75      0.81      1.00;
+      0.05      0.80      1.00;
+     -2.18      0.82      1.00;
+      3.88      0.81      1.00;
+      0.73      0.77      1.00;
+      5.70      0.77      1.00;
+     -0.69      0.83      1.00;
+     -8.14      0.87      1.00;
+      5.37      0.99      1.00;
+      1.87      0.95      1.00;
+      5.76      0.80      1.00;
+     -0.79      0.89      1.00;
+    -13.23      1.21      1.00;
+      3.97      1.26      1.00;
+      5.20      0.81      1.00;
+      3.16      0.61      1.00;
+      6.41      0.53      1.00;
+      1.72      0.64      1.00;
+      2.20      0.75      1.00;
+      1.06      0.95      1.00;
+      9.54      0.96      1.00;
+     -4.75      1.31      1.00;
+     -5.05      1.04      1.00;
+      0.48      1.07      1.00;
+      3.41      1.21      1.00;
+     -2.21      1.08      1.00;
+      0.21      1.15      1.00;
+     -2.37      1.35      1.00;
+     -1.55      1.24      1.00;
+     -6.91      1.28      1.00;
+     -7.62      1.24      1.00;
+      4.81      1.21      1.00;
+      3.51      1.07      1.00;
+     -3.68      0.87      1.00;
+     -3.42      0.80      1.00;
+     -6.03      0.92      1.00;
+     -1.99      0.98      1.00;
+      3.20      1.13      1.00;
+     -3.88      1.06      1.00;
+     -3.35      0.96      1.00;
+     -3.10      1.05      1.00;
+     11.14      0.76      1.00;
+      1.17      0.51      1.00;
+     11.27      0.59      1.00;
+      4.56      0.63      1.00;
+      0.78      0.67      1.00;
+      3.50      0.69      1.00;
+      2.40      0.62      1.00;
+      2.84      0.63      1.00;
+      6.71      0.71      1.00;
+      0.63      0.69      1.00;
+      3.11      0.67      1.00;
+     -3.90      0.74      1.00;
+     -0.41      0.76      1.00;
+      0.85      0.76      1.00;
+     -3.56      0.76      1.00;
+      2.26      0.70      1.00;
+     -1.78      0.73      1.00;
+     -2.06      0.76      1.00;
+     -4.62      0.71      1.00;
+      0.61      0.73      1.00;
+     -0.56      0.81      1.00;
+     -6.01      0.78      1.00;
+      1.59      0.75      1.00;
+     -2.88      0.82      1.00;
+     10.44      0.83      1.00;
+     -0.82      0.86      1.00;
+     -1.01      1.00      1.00;
+     -1.80      0.73      1.00;
+      1.73      0.64      1.00;
+      7.92      0.65      1.00;
+      1.11      0.58      1.00;
+     -0.79      0.62      1.00;
+     -0.94      0.72      1.00;
+      4.92      0.66      1.00;
+      1.16      0.55      1.00;
+     -0.65      0.62      1.00;
+     -1.03      0.55      1.00;
+     -4.58      0.60      1.00;
+      3.79      0.65      1.00;
+      6.31      0.61      1.00;
+      3.66      0.65      1.00;
+      0.42      0.56      1.00;
+      6.72      0.53      1.00;
+      4.79      0.60      1.00;
+     -1.31      0.52      1.00;
+      4.59      0.49      1.00;
+      0.90      0.52      1.00;
+     -6.49      0.52      1.00;
+      6.16      0.46      1.00;
+     -8.35      0.45      1.00;
+      4.47      0.46      1.00;
+      1.12      0.39      1.00;
+     -3.13      0.49      1.00;
+     12.43      0.42      1.00;
+      4.36      0.43      1.00;
+      1.90      0.47      1.00;
+     -2.14      0.44      1.00 ]
+
+e = [ 5.19346;
+      2.22296;
+      3.82622;
+      2.93541;
+      0.70453;
+      1.98955;
+      0.56655;
+      1.54306;
+      0.57363;
+     -0.13346;
+      1.68565;
+      9.24558;
+      7.39317;
+      2.45484;
+     -4.83724;
+     -0.74742;
+     -0.22633;
+      0.51537;
+      2.82822;
+      9.35943;
+      5.38263;
+      6.59691;
+      1.79080;
+      0.54506;
+      4.06661;
+     -4.18808;
+      1.69414;
+      5.40963;
+      2.92762;
+     -2.67351;
+     -2.74613;
+     -3.52752;
+     -4.27775;
+      2.88245;
+     -8.12566;
+      0.08760;
+      2.20495;
+     -0.45509;
+     -0.06959;
+      1.51330;
+      1.37878;
+     -2.33046;
+      0.35730;
+    -11.41189;
+     -0.48083;
+      1.08917;
+      2.52333;
+      1.73218;
+     15.60426;
+     -0.15961;
+     -0.22174;
+     -1.19174;
+     13.92022;
+      1.16044;
+      0.08087;
+     -5.81547;
+     -5.14199;
+     -5.69918;
+     -0.36585;
+     -3.70294;
+      2.83687;
+     -1.37631;
+     -2.55403;
+     -4.10796;
+      1.32092;
+     -1.64056;
+     -4.09963;
+     -4.85728;
+     -1.45496;
+     -3.12612;
+     -4.36648;
+     -3.23379;
+      5.22893;
+      5.06730;
+     -0.63041;
+     -1.70377;
+     -5.71335;
+     -2.76831;
+      2.20443;
+     -0.10174;
+     -2.36727;
+     -5.49859;
+     -4.78093;
+     -3.36665;
+      5.12027;
+     -3.03180;
+     -2.52035;
+      3.10913;
+     -3.51388;
+     -0.70389;
+     -1.01963;
+     -8.78989;
+      2.95170;
+     -3.95464;
+     -4.01258;
+     -1.92315;
+     -4.94209;
+      5.29138;
+      1.28435;
+      0.28832 ]
+
+    atol = 1e-4
+    w_test = WhiteTest(X, e, type = :linear)
+    @test pvalue(w_test) ≈ 0.1287 atol=atol
+
+    w_test = WhiteTest(X, e, type = :linear_and_squares)
+    @test pvalue(w_test) ≈ 0.2774 atol=atol
+
+    w_test = WhiteTest(X, e)
+    w_pval = pvalue(w_test)
+    @test w_pval ≈ 0.3458 atol=atol
+
+    refstr = """
+    White's (or Breusch-Pagan's) test for heteroskedasticity
+    --------------------------------------------------------
+    Population details:
+        parameter of interest:   T*R2
+        value under h_0:         0
+        point estimate:          $(w_test.lm)
+
+    Test summary:
+        outcome with 95% confidence: fail to reject h_0
+        one-sided p-value:           $(round(w_pval,digits=4))
+
+    Details:
+        T*R^2 statistic:        $(w_test.lm)
+        degrees of freedom:     $(w_test.dof)
+        type:                   $(w_test.type)
+    """
+
+    buffer = IOBuffer()
+    show(buffer, w_test)
+    str = String(take!(buffer))
+    @test str == refstr
+
+    @test_throws DimensionMismatch WhiteTest([rand(3) ones(3)], rand(4))
+
+    @test_throws ArgumentError WhiteTest(ones(4,1), rand(4))
+
+    @test_throws ArgumentError WhiteTest([zeros(4) rand(4)], rand(4))
+
+    @test_throws ArgumentError WhiteTest(rand(4,2), rand(4))
+
+    bp_test = BreuschPaganTest(X, e)
+    @test pvalue(bp_test) ≈ 0.1287 atol=atol
+    show(IOBuffer(), bp_test)
+end