change tolerance to 1e-6 (#166)

* Update fit.jl

* Update fit.jl

* Update Project.toml

* Create ci.yml

* Update README.md

* update

* update README

* Update README.md

* rework on Readme

Author: matthieugomez

.github/workflows/ci.yml

@@ -0,0 +1,46 @@
+name: CI
+on:
+  pull_request:
+    branches:
+      - master
+  push:
+    branches:
+      - master
+    tags: '*'
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.3' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
+          - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
+

.gitignore

@@ -0,0 +1,2 @@
+
+benchmark/.sublime2Terminal.jl

Project.toml

@@ -1,6 +1,6 @@
 name = "FixedEffectModels"
 uuid = "9d5cd8c9-2029-5cab-9928-427838db53e3"
-version = "1.5.1"
+version = "1.6.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

README.md

@@ -1,4 +1,4 @@
-[![Build Status](https://travis-ci.com/FixedEffects/FixedEffectModels.jl.svg?branch=master)](https://travis-ci.com/FixedEffects/FixedEffectModels.jl)
+[![Build status](https://github.com/FixedEffects/FixedEffectModels.jl/workflows/CI/badge.svg)](https://github.com/FixedEffects/FixedEffectModels.jl/actions)
 
 This package estimates linear models with high dimensional categorical variables and/or instrumental variables.
 
@@ -11,7 +11,7 @@ The objective of the package is similar to the Stata command [`reghdfe`](https:/
 ![benchmark](http://www.matthieugomez.com/files/fixedeffectmodels_benchmark.png)
 
 
-Performances are roughly similar to the newer R function [`feols`](https://cran.r-project.org/web/packages/fixest/fixest.pdf) (note: use `tol = 1e-6, drop_singletons = false` to match the default options of `feols`). The main difference is that `FixedEffectModels` can also run the demeaning operation on a GPU (with `method = :gpu`).
+Performances are roughly similar to the newer R function [`feols`](https://cran.r-project.org/web/packages/fixest/fixest.pdf). The main difference is that `FixedEffectModels` can also run the demeaning operation on a GPU (with `method = :gpu`).
 
 ## Syntax
 
@@ -55,23 +55,29 @@ reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)), Vcov.cluster(:State), weig
 	reg(df, term(:Sales) ~ term(:NDI) + fe(:State) + fe(:Year))
 	```
 
-- Standard errors are indicated with the prefix `Vcov` (with the package [Vcov](http://github.com/matthieugomez/Vcov.jl))
+- The option `contrasts` specifies that a column should be understood as a set of dummy variables:
 	```julia
-	Vcov.robust()
-	Vcov.cluster(:State)
-	Vcov.cluster(:State, :Year)
+	reg(df, @formula(Sales ~ Price + Year); contrasts = Dict(:Year => DummyCoding()))
 	```
+	You can specify different base levels 
+	```julia
+	reg(df, @formula(Sales ~ Price + Year); contrasts = Dict(:Year => DummyCoding(base = 80)))
+	```
+
+
 - The option `weights` specifies a variable for weights
 	```julia
 	weights = :Pop
 	```
-- The option `contrasts` specifies particular contrasts for a dummy variable in the formula, e.g.
+
+- Standard errors are indicated with the prefix `Vcov` (with the package [Vcov](http://github.com/matthieugomez/Vcov.jl))
 	```julia
-	reg(df, @formula(Sales ~ Year); contrasts = Dict(:Year => DummyCoding(base = 80)))
+	Vcov.robust()
+	Vcov.cluster(:State)
+	Vcov.cluster(:State, :Year)
 	```
-- The option `save` can be set to one of the following:  `none` (default) to save nothing `:residuals` to save residuals, `:fe` to save fixed effects. You can access the result with `residuals()` and `fe()`
 
--
+- The option `save` can be set to one of the following:  `none` (default) to save nothing `:residuals` to save residuals, `:fe` to save fixed effects. You can access the result with `residuals()` and `fe()`
 
 - The option `method` can be set to one of the following: `:cpu`, `:gpu` (see Performances below).
 
@@ -93,10 +99,11 @@ You may use [RegressionTables.jl](https://github.com/jmboehm/RegressionTables.jl
 
 ## Performances
 
-### MultiThreads
-By default, `FixedEffectModels` uses as many threads as `Threads.nthreads()`.  Use the option `nthreads` to select the number of threads to use in the estimation. Default to `Threads.nthreads()`.
 
+### MultiThreads
+`FixedEffectModels` is multi-threaded. Use the option `nthreads` to select the number of threads to use in the estimation (defaults to `Threads.nthreads()`). That being said, multithreading does not usually make a big difference.
 
+### GPU
 The package has support for GPUs (Nvidia) (thanks to Paul Schrimpf). This can make the package an order of magnitude faster for complicated problems.
 
 To use GPU, run `using CUDA` before `using FixedEffectModels`. Then, estimate a model with `method = :gpu`. For maximum speed, set the floating point precision to `Float32` with `double_precision = false`.

benchmark/benchmark.jl

@@ -1,4 +1,4 @@
-using DataFrames, FixedEffectModels, Random
+using DataFrames, FixedEffectModels, Random, CategoricalArrays
 
 # Very simple setup
 N = 10000000
@@ -17,9 +17,9 @@ df = DataFrame(id1 = id1, id2 = id2, x1 = x1, x2 = x2, y = y)
 @time reg(df, @formula(y ~ x1 + x2),  Vcov.cluster(:id2))
 # 0.759032 seconds (631 allocations: 728.437 MiB, 4.89% gc time)
 @time reg(df, @formula(y ~ x1 + x2 + fe(id1)))
-# 1.452800 seconds (231.14 k allocations: 908.415 MiB, 3.40% gc time, 9.72% compilation time)
+# 1.316560 seconds (230.66 k allocations: 908.386 MiB, 4.28% gc time, 11.84% compilation time)
 @time reg(df, @formula(y ~ x1 + x2 + fe(id1)), Vcov.cluster(:id1))
-# 1.785358 seconds (231.43 k allocations: 952.060 MiB, 8.22% compilation time)
+# 1.501165 seconds (230.94 k allocations: 952.029 MiB, 0.84% gc time, 10.66% compilation time)
 @time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
 # 3.058090 seconds (331.56 k allocations: 1.005 GiB, 2.08% gc time, 5.89% compilation time)
 
@@ -35,6 +35,22 @@ x2 =  cos.(id1) +  sin.(id2) + randn(N)
 y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
 df = DataFrame(id1 = id1, id2 = id2, x1 = x1, x2 = x2, y = y)
 @time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))
-# 15.738214 seconds (929.97 k allocations: 147.737 MiB, 0.95% gc time)
+#   3.190288 seconds (393.89 k allocations: 109.614 MiB, 3.63% gc time, 9.75% compilation time)
 
 
++# fixest
+n = 10_000_000
+nb_dum = [div(n,20), floor(Int, sqrt(n)), floor(Int, n^.33)]
+N = nb_dum.^3
+id1 = categorical(rand(1:nb_dum[1], n))
+id2 = categorical(rand(1:nb_dum[2], n))
+id3 = categorical(rand(1:nb_dum[3], n))
+X1 = rand(n)
+ln_y = 3 .* X1 .+ rand(n) 
+df = DataFrame(X1 = X1, ln_y = ln_y, id1 = id1, id2 = id2, id3 = id3)
+@time reg(df, @formula(ln_y ~ X1 + fe(id1)), Vcov.cluster(:id1))
+# 1.058386 seconds (219.24 k allocations: 751.826 MiB, 13.05% compilation time)
+@time reg(df, @formula(ln_y ~ X1 + fe(id1) + fe(id2)), Vcov.cluster(:id1))
+# 2.895223 seconds (284.28 k allocations: 870.782 MiB, 1.73% gc time, 5.97% compilation time)
+@time reg(df, @formula(ln_y ~ X1 + fe(id1) + fe(id2) + fe(id3)), Vcov.cluster(:id1))
+# 3.191483 seconds (380.10 k allocations: 991.134 MiB, 6.88% compilation time)

src/FixedEffectModel.jl

@@ -191,7 +191,6 @@ function top(x::FixedEffectModel)
         out = vcat(out, 
             ["R2 within" format_scientific(x.r2_within);
            "Iterations" sprint(show, x.iterations, context = :compact => true);
-             "Converged" sprint(show, x.converged, context = :compact => true);
              ])
     end
     return out

src/fit.jl

@@ -13,7 +13,7 @@ Estimate a linear model with high dimensional categorical variables / instrument
 * `method::Symbol`: A symbol for the method. Default is :cpu. Alternatively,  :gpu requires `CuArrays`. In this case, use the option `double_precision = false` to use `Float32`.
 * `nthreads::Integer` Number of threads to use in the estimation. If `method = :cpu`, defaults to `Threads.nthreads()`. If `method = :gpu`, defaults to 256.
 * `double_precision::Bool`: Should the demeaning operation use Float64 rather than Float32? Default to true.
-* `tol::Real` Tolerance. Default to 1e-8 if `double_precision = true`, 1e-6 otherwise.
+* `tol::Real` Tolerance. Default to 1e-6.
 * `maxiter::Integer = 10000`: Maximum number of iterations
 * `drop_singletons::Bool = true`: Should singletons be dropped?
 * `progress_bar::Bool = true`: Should the regression show a progressbar
@@ -50,7 +50,7 @@ function reg(
     @nospecialize(method::Symbol = :cpu),
     @nospecialize(nthreads::Integer = method == :cpu ? Threads.nthreads() : 256),
     @nospecialize(double_precision::Bool = true),
-    @nospecialize(tol::Real = double_precision ? 1e-8 : 1e-6),
+    @nospecialize(tol::Real = 1e-6),
     @nospecialize(maxiter::Integer = 10000),
     @nospecialize(drop_singletons::Bool = true),
     @nospecialize(progress_bar::Bool = true),
@@ -447,5 +447,6 @@ function reg(
     if esample == Colon()
         esample = trues(N)
     end
+
     return FixedEffectModel(coef, matrix_vcov, vcov, nclusters, esample, residuals2, augmentdf, fekeys, coef_names, response_name, formula_origin, formula, contrasts, nobs, dof_residual_,  rss, tss_total, r2, adjr2, F, p, iterations, converged, r2_within, F_kp, p_kp)
 end

test/fit.jl

@@ -79,6 +79,7 @@ x = reg(df, m)
 m = @formula y ~ z1 + x1&x2&fe(id1)
 x = reg(df, m)
 @test coef(x)  ≈  [122.98713, 0.30933] atol = 1e-4
+# SSR does not work well here
 m = @formula y ~ z1 + (x1&x2)*fe(id1)
 x = reg(df, m)
 @test coef(x) ≈   [0.421406, 0.0] atol = 1e-4
@@ -612,11 +613,13 @@ x = reg(df, m, weights = :w)
 @test coef(x) ≈  [- 0.043475472188120416] atol = 1e-3
 
 ## the last two ones test an ill conditioned model matrix
+# SSR does not work well here
 m = @formula y ~ x1 + fe(id1) + fe(id1)&id2
 x = reg(df, m)
 @test coef(x)  ≈   [- 0.122354] atol = 1e-4
 @test x.iterations <= 30
 
+# SSR does not work well here
 m = @formula y ~ x1 + fe(id1) + fe(id1)&id2
 x = reg(df, m, weights = :w)
 @test coef(x) ≈ [- 0.11752306001586807] atol = 1e-4