Repository: JuliaStats/HypothesisTests.jl
Pull request #167 (2019-10-28 19:55:37)

Adding homoscedasticity test
Hi, I'm interested in contributing with a homoscedasticity test. The main reference is Durbin & Koopman, Time Series Analysis by State Space Methods, 2nd ed., chapter 2.12.1.

This is my first time using the package so please tell me if something is off. I tried to make it as similar as possible to how the other tests were implemented.

raphaelsaavedra: Hi @nalimilan,

Thanks for the review. Unfortunately, I don't have a proper name for the test. My reference is Durbing & Koopman and they only mention it as H test (name they gave to the test statistic). 

<img width="974" alt="Screen Shot 2019-11-02 at 16 33 58" src="https://user-images.githubusercontent.com/29098169/68076042-afe95a80-fd8e-11e9-8d41-4570a11bee32.png">

I will try to find out if the test has a proper name.
nalimilan: It sounds weird that they would propose a new test with so little developments. It would be good to find a longer reference and a standard name.
raphaelsaavedra: Well it's mostly because it's simply an F-test that is meant to be applied to residuals. If you find it acceptable, I can rename it to `HomoscedasticityFTest` or something like that.
mschauer: It would be good to provide the F test independently from the application to time series and make the HomoscedasticityFTest test a thin wrapper. (I hope that is not too much feature creep) 
raphaelsaavedra: Very good point @mschauer. I am changing the struct to `FTest` so that we have a more general form that can be used to test the variances of two different samples and then `HomoscedasticityFTest` will simply be a particular case of `FTest`.
nalimilan: There are several kinds of F-tests, as can be seen on Wikipedia. Better use a more specific name than `FTest` to avoid confusion with e.g. `ftest` in GLM. We currently have `EqualVarianceTTest`, but unfortunately "equal variance" there refers to the assumption of the t-test, not to the hypothesis which is tested. So maybe use just `VarianceFTest`? `EqualityOfVariancesFTest` would be more explicit but quite verbose.

Anyway, probably better list this test under "Parametric tests", and put only the homoscedasticity test under time series.

Regarding `HomoscedasticityFTest`, I don't think it's appropriate to have a function which looks like a constructor, but returns an object of a different type. It should either use a separate type, or make it a constructor of `VarianceFTest`. Not sure which approach is better. I guess the advantage of the latter is that it would allow printing the value of the `h` argument, and it would be simpler for the docs (list it with time series). The drawback is of course the increased complexity, but the code is quite short. What do you think?
raphaelsaavedra: I am making the requested changes. I agree with your point regarding `HomoscedasticityFTest`, but if we change it to a constructor of `VarianceFTest` I'm not sure how to present the different docstrings for methods of the same name under parametric and time series. Do you know how to do this? 
nalimilan: You can give the signature of each method instead of just its name. But I'm still not completely sure whether it's better to use a single type or two different types...
raphaelsaavedra: I think having a single constructor for `VarianceFTest` is pretty intuitive. One method that receives two vectors and another one that receives one vector and the number of observations to consider in each "side".
nalimilan: OK. Though I guess both could be listed under "Parametric tests", as it could apply to series other than time series.
raphaelsaavedra: See if it's ok now :) don't hesitate to ping me to make additional changes
raphaelsaavedra: I implemented your suggestions @nalimilan, please tell me if it's okay now.
mschauer: I understand that this is opposite to what you wanted to do in the first place, but I would be most comfortable with just adding the `VarianceFTest(y1::AbstractVector{<: Real}, y2::AbstractVector{<: Real})` and defer any time series applications to that to user space or more specialised packages as
```
    y1 = @view y[n-h+1:n]
    y2 = @view y[1:h]
```
is easily done on the  callers site and that `VarianceFTest` with one argument means splitting the vector is a bit too magic. As a compromise, having a *method* `homoscedasticityFTest(y::AbstractVector{<: Real}; firstlast::Integer = div(length(y), 2)) = ...` returning a VarianceFTest would be okay with me.
raphaelsaavedra: Thinking about it, I have to agree with @mschauer. The homoscedasticity test is such a simple case of `VarianceFTest` that it does not justify its own constructor. I don't even think having its own method is necessary. If @nalimilan agrees, I can just remove the homoscedasticity version of `VarianceFTest` and anyone who wants to run that test can simply do as mentioned by @mschauer.
nalimilan: Fine with me.

Patch:
diff --git a/docs/src/parametric.md b/docs/src/parametric.md
index 10f1452..2e463b2 100644
--- a/docs/src/parametric.md
+++ b/docs/src/parametric.md
@@ -28,3 +28,8 @@ OneSampleZTest
 EqualVarianceZTest
 UnequalVarianceZTest
 ```
+
+## F-test
+```@docs
+VarianceFTest
+```
\ No newline at end of file
diff --git a/docs/src/time_series.md b/docs/src/time_series.md
index ed10253..bfe77d7 100644
--- a/docs/src/time_series.md
+++ b/docs/src/time_series.md
@@ -22,4 +22,4 @@ JarqueBeraTest
 ## Augmented Dickey-Fuller test
 ```@docs
 ADFTest
-```
+```
\ No newline at end of file
diff --git a/src/f.jl b/src/f.jl
new file mode 100644
index 0000000..2e1ed06
--- /dev/null
+++ b/src/f.jl
@@ -0,0 +1,77 @@
+# f.jl
+# F-tests
+#
+# Copyright (C) 2019   Raphael Saavedra
+#
+# 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 VarianceFTest
+
+struct VarianceFTest <: HypothesisTest
+    n_x::Int   # size of sample 1
+    n_y::Int   # size of sample 2
+    df_x::Int  # degrees of freedom 1
+    df_y::Int  # degrees of freedom 2
+    F::Real    # test statistic
+end
+
+"""
+    VarianceFTest(y1::AbstractVector{<: Real}, y2::AbstractVector{<: Real})
+
+Perform an F-test of the null hypothesis that two real-valued vectors `y1` and `y2` have equal variances.
+
+Implements: [`pvalue`](@ref)
+
+# References
+
+  * George E. P. Box, "Non-Normality and Tests on Variances", Biometrika 40 (3/4): 318–335, 1953.
+
+# External links
+
+  * [F-test of equality of variances on Wikipedia](https://en.wikipedia.org/wiki/F-test_of_equality_of_variances)
+"""
+function VarianceFTest(y1::AbstractVector{<: Real}, y2::AbstractVector{<: Real})
+    n1, n2 = length(y1), length(y2)
+    F = var(y1) / var(y2)
+    return VarianceFTest(n1, n2, n1-1, n2-1, F)
+end
+
+testname(::VarianceFTest) = "Variance F-test"
+population_param_of_interest(x::VarianceFTest) = ("variance ratio", 1.0, x.F)
+default_tail(test::VarianceFTest) = :both
+
+function show_params(io::IO, x::VarianceFTest, ident)
+    println(io, ident, "number of observations: [$(x.n_x), $(x.n_y)]")
+    println(io, ident, "F statistic:            $(x.F)")
+    println(io, ident, "degrees of freedom:     [$(x.df_x), $(x.df_y)]")
+end
+
+function pvalue(x::VarianceFTest; tail=:both)
+    dist = FDist(x.df_x, x.df_y)
+    if tail == :both
+        return 1 - 2*abs(cdf(dist, x.F) - 0.5)
+    elseif tail == :right
+        return 1 - cdf(dist, x.F)
+    elseif tail == :left
+        return cdf(dist, x.F)
+    else
+        throw(ArgumentError("tail=$(tail) is invalid"))
+    end
+end
\ No newline at end of file


Test patch:
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index 0963fc0..12b45a1 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -188,4 +188,5 @@ include("permutation.jl")
 include("hotelling.jl")
 include("bartlett.jl")
 include("wald_wolfowitz.jl")
+include("f.jl")
 end
diff --git a/test/f.jl b/test/f.jl
new file mode 100644
index 0000000..a0ab58f
--- /dev/null
+++ b/test/f.jl
@@ -0,0 +1,49 @@
+using HypothesisTests, Test
+using HypothesisTests: default_tail
+
+@testset "F-tests" begin
+    @testset "Basic variance F-test" begin
+        Random.seed!(12)
+        y1_h0 = 4 .+ randn(500)
+        y2_h0 = 4 .+ randn(400)
+
+        t = VarianceFTest(y1_h0, y2_h0)
+
+        @test t.n_x == 500
+        @test t.n_y == 400
+        @test t.df_x == 499
+        @test t.df_y == 399
+        @test t.F ≈ 0.974693 rtol=1e-5
+        @test pvalue(t) ≈ 0.784563 rtol=1e-5
+        @test pvalue(t, tail=:left) ≈ 0.392281 rtol=1e-5
+        @test pvalue(t, tail=:right) ≈ 0.607718 rtol=1e-5
+        @test default_tail(t) == :both
+
+        t = VarianceFTest(y2_h0, y1_h0)
+
+        @test t.n_x == 400
+        @test t.n_y == 500
+        @test t.df_x == 399
+        @test t.df_y == 499
+        @test t.F ≈ 1.025963 rtol=1e-5
+        @test pvalue(t) ≈ 0.784563 rtol=1e-5
+        @test pvalue(t, tail=:right) ≈ 0.392281 rtol=1e-5
+        @test pvalue(t, tail=:left) ≈ 0.607718 rtol=1e-5
+        @test default_tail(t) == :both
+
+        y1_h1 = 0.7*randn(200)
+        y2_h1 = 1.3*randn(120)
+
+        t = VarianceFTest(y1_h1, y2_h1)
+
+        @test t.n_x == 200
+        @test t.n_y == 120
+        @test t.df_x == 199
+        @test t.df_y == 119
+        @test t.F ≈ 0.367547 rtol=1e-5
+        @test pvalue(t) < 1e-8
+        @test default_tail(t) == :both
+        @test pvalue(t, tail=:left) < 1e-8
+        @test pvalue(t, tail=:right) > 1.0 - 1e-8
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7b2fe5c..ade795d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -34,3 +34,4 @@ include("t.jl")
 include("wald_wolfowitz.jl")
 include("wilcoxon.jl")
 include("z.jl")
+include("f.jl")
\ No newline at end of file