Refractor of struct

src/Econometrics.jl

@@ -16,7 +16,7 @@ module Econometrics
                          Diagonal, Hermitian, I, LowerTriangular, qr,
                          UpperTriangular
     using Optim: hessian!, optimize, minimizer, TwiceDifferentiable
-    using Parameters: @unpack
+    using Parameters: @unpack, @pack!
     using Printf: @sprintf
     using StatsBase: aic, aicc, bic, harmmean, FrequencyWeights, CoefTable,
                      ConvergenceException, Weights
@@ -35,13 +35,14 @@ module Econometrics
                       predict!, dof_residual, RegressionModel, params
     import StatsModels: hasintercept, implicit_intercept
     foreach(file -> include(joinpath(dirname(@__DIR__), "src", "$file.jl")),
-            ["structs", "transformations", "formula", "solvers", "main", "statsbase", "wald"])
+            ["structs", "transformations", "formula", "main", "solvers", "statsbase", "wald"])
     export @formula, DummyCoding, aic, aicc, bic, coef, coefnames,
            coeftable, confint, deviance, islinear, nulldeviance, loglikelihood,
            nullloglikelihood, score, nobs, dof, mss, rss, informationmatrix, vcov, stderror,
            weights, isfitted, fit, fit!, r2, adjr2, fitted, response, meanresponse,
            modelmatrix, leverage, residuals, predict, predict!, dof_residual,
            params,
            EconometricModel, absorb, BetweenEstimator, RandomEffectsEstimator,
+           ContinuousResponse,
            OIM, HC0, HC1, HC2, HC3, HC4
 end

src/formula.jl

@@ -127,7 +127,7 @@ function decompose(f::FormulaTerm,
     end
     if !isa(iv.lhs, InterceptTerm)
         z, Z = modelcols(iv, data)
-        z = isa(z, Tuple) ? hcat(z...) : z
+        z = isa(z, Tuple) ? hcat(z...) : convert(Matrix{Float64}, reshape(z, length(z), 1))
     else
         z = zeros(0, 0)
         Z = zeros(0, 0)

src/main.jl

@@ -27,15 +27,48 @@ mutable struct EconometricModel{E<:ModelEstimator,
     data::T
     X::Matrix{Float64}
     y::Y
-    w::W
+    wts::W
     β::Vector{Float64}
     Ψ::Hermitian{Float64,Matrix{Float64}}
     ŷ::Ŷ
+    z::Matrix{Float64}
+    Z::Matrix{Float64}
     vars::N
     iv::Int
     vce::VC
+    function EconometricModel(estimator::Type{<:Union{EconometricModel,ModelEstimator}},
+                              f::FormulaTerm,
+                              data;
+                              contrasts::Dict{Symbol} = Dict{Symbol,Union{<:AbstractContrasts,<:AbstractTerm}}(),
+                              wts::Union{Nothing,Symbol} = nothing,
+                              panel::Union{Nothing,Symbol} = nothing,
+                              time::Union{Nothing,Symbol} = nothing,
+                              vce::VCE = OIM)
+        data, exogenous, iv, estimator, X, y, z, Z, wts =
+            decompose(deepcopy(f), data, contrasts, wts, panel, time, estimator, vce)
+        if isa(estimator, Union{NominalResponse, OrdinalResponse})
+            @unpack categories = estimator
+            y = [ findfirst(isequal(x), categories) for x ∈ y ]
+            @assert length(categories) > 2
+        end
+        if isa(estimator, NominalResponse)
+            ŷ = zeros(0, 0)
+        else
+            ŷ = zeros(0)
+        end
+        wts = FrequencyWeights(collect(wts))
+        vars = (coefnames(exogenous.lhs),
+                convert(Vector{String}, vcat(coefnames(exogenous.rhs), coefnames(iv.lhs))))
+        new{typeof(estimator), typeof(f), typeof(data), typeof(y), typeof(wts), typeof(ŷ), typeof(vars), typeof(vce)}(
+            estimator, f, data, X, y, wts, zeros(0), Hermitian(zeros(0, 0)), ŷ,
+            z, Z, vars, size(Z, 2), vce)
+    end
 end
 function show(io::IO, obj::EconometricModel{<:LinearModelEstimators})
+    if !isfitted(obj)
+        println(io, "Model has not been fitted.")
+        return obj
+    end
     show(io, obj.estimator)
     println(io, @sprintf("Number of observations: %i", nobs(obj)))
     println(io, @sprintf("Null Loglikelihood: %.2f", nullloglikelihood(obj)))
@@ -50,6 +83,10 @@ function show(io::IO, obj::EconometricModel{<:LinearModelEstimators})
     show(io, coeftable(obj))
 end
 function show(io::IO, obj::EconometricModel{<:ContinuousResponse})
+    if !isfitted(obj)
+        println(io, "Model has not been fitted.")
+        return obj
+    end
     show(io, obj.estimator)
     ℓℓ₀ = nullloglikelihood(obj)
     ℓℓ = loglikelihood(obj)
@@ -78,6 +115,10 @@ function show(io::IO, obj::EconometricModel{<:ContinuousResponse})
     show(io, coeftable(obj))
 end
 function show(io::IO, obj::EconometricModel{<:NominalResponse})
+    if !isfitted(obj)
+        println(io, "Model has not been fitted.")
+        return obj
+    end
     show(io, obj.estimator)
     ℓℓ₀ = nullloglikelihood(obj)
     ℓℓ = loglikelihood(obj)
@@ -97,6 +138,10 @@ function show(io::IO, obj::EconometricModel{<:NominalResponse})
     show(io, coeftable(obj))
 end
 function show(io::IO, obj::EconometricModel{<:OrdinalResponse})
+    if !isfitted(obj)
+        println(io, "Model has not been fitted.")
+        return obj
+    end
     show(io, obj.estimator)
     ℓℓ₀ = nullloglikelihood(obj)
     ℓℓ = loglikelihood(obj)
@@ -118,15 +163,14 @@ end
 function fit(estimator::Type{<:Union{EconometricModel,ModelEstimator}},
              f::FormulaTerm,
              data;
-             contrasts::Dict{Symbol} = Dict{Symbol,Union{<:AbstractContrasts,<:AbstractTerm}}(),
-             weights::Union{Nothing,Symbol} = nothing,
-             panel::Union{Nothing,Symbol} = nothing,
-             time::Union{Nothing,Symbol} = nothing,
-             vce::VCE = OIM)
-    data, exogenous, iv, estimator, X, y, z, Z, wts =
-        decompose(deepcopy(f), data, contrasts, weights, panel, time, estimator, vce)
-    X, y, β, Ψ, ŷ, wts, piv = solve(estimator, X, y, z, Z, wts)
-    vars = (coefnames(exogenous.lhs),
-            convert(Vector{String}, vcat(coefnames(exogenous.rhs), coefnames(iv.lhs))[piv]))
-    EconometricModel(estimator, f, data, X, y, wts, β, Ψ, ŷ, vars, size(Z, 2), vce)
+             fit = true,
+             kw...)
+    model = EconometricModel(estimator, f, data;
+                             contrasts = get(kw, :contrasts, Dict{Symbol,Union{<:AbstractContrasts,<:AbstractTerm}}()),
+                             wts = get(kw, :wts, nothing),
+                             panel = get(kw, :panel, nothing),
+                             time = get(kw, :time, nothing),
+                             vce = get(kw, :vce, OIM))
+    fit && fit!(model)
+    model
 end

src/solvers.jl

@@ -15,6 +15,50 @@ For example, Swamy-Arora random effects model for longitudinal data.
                       z::AbstractVecOrMat{<:Number},
                       Z::AbstractMatrix{<:Number},
                       wts::AbstractVector)
+	if !isa(wts, Weights)
+        wts = FrequencyWeights(wts)
+    else
+        isa(wts, FrequencyWeights) || throw(ArgumentError("Only frequency weights are supported."))
+    end
+    X = transform(estimator, X, wts)
+    y = transform(estimator, y, wts)
+    z = transform(estimator, z, wts)
+    Z = transform(estimator, Z, wts)
+    wts = transform(estimator, FrequencyWeights(wts))
+    if !isempty(z)
+        Z̃ = hcat(X, Z)
+        F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true, check = false)
+        bkr = count(x -> abs(x) ≥ √eps(), F.D)
+		if bkr < size(F, 2)
+            lin_ind = sort!(invperm(F.p)[1:bkr])
+            count(x -> lin_ind > size(X, 2) + 1) ≥ size(z, 2) ||
+                throw(ArgumentError("Insufficient number of instruments."))
+            Z̃ = Z̃[:,lin_ind]
+            F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true)
+        else
+            lin_ind = collect(1:size(F, 2))
+        end
+        γ = F \ (Z̃' * Diagonal(wts) * z)
+        X̃ = hcat(X, Z̃ * γ)
+    else
+        X̃ = X
+    end
+    F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true, check = false)
+    bkr = count(x -> abs(x) ≥ √eps(), F.D)
+    if bkr < size(F, 2)
+        lin_ind = sort!(invperm(F.p)[1:bkr])
+        X̃ = convert(Matrix{Float64}, X̃[:,lin_ind])
+        F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true)
+    else
+        lin_ind = collect(1:size(F, 2))
+    end
+    β = F \ (X̃' * Diagonal(wts) * y)
+    Ψ = Hermitian(inv(F))
+    ŷ = isempty(z) ? X̃ * β : hcat(X, z)[:,lin_ind] * β
+	X̃, y, β, Ψ, ŷ, wts, lin_ind
+end
+@views function fit!(obj::EconometricModel{<:LinearModelEstimators})
+	@unpack estimator, X, y, z, Z, wts = obj
     if !isa(wts, Weights)
         wts = FrequencyWeights(wts)
     else
@@ -24,10 +68,10 @@ For example, Swamy-Arora random effects model for longitudinal data.
     y = transform(estimator, y, wts)
     z = transform(estimator, z, wts)
     Z = transform(estimator, Z, wts)
-    w = transform(estimator, FrequencyWeights(wts))
+    wts = transform(estimator, FrequencyWeights(wts))
     if !isempty(z)
         Z̃ = hcat(X, Z)
-        F = bunchkaufman(Hermitian(Z̃' * Diagonal(w) * Z̃), true, check = false)
+        F = bunchkaufman(Hermitian(Z̃' * Diagonal(wts) * Z̃), true, check = false)
         bkr = count(x -> abs(x) ≥ √eps(), F.D)
 		if bkr < size(F, 2)
             lin_ind = sort!(invperm(F.p)[1:bkr])
@@ -38,41 +82,31 @@ For example, Swamy-Arora random effects model for longitudinal data.
         else
             lin_ind = collect(1:size(F, 2))
         end
-        γ = F \ (Z̃' * Diagonal(w) * z)
+        γ = F \ (Z̃' * Diagonal(wts) * z)
         X̃ = hcat(X, Z̃ * γ)
     else
         X̃ = X
     end
-    F = bunchkaufman(Hermitian(X̃' * Diagonal(w) * X̃), true, check = false)
+    F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true, check = false)
     bkr = count(x -> abs(x) ≥ √eps(), F.D)
     if bkr < size(F, 2)
         lin_ind = sort!(invperm(F.p)[1:bkr])
         X̃ = convert(Matrix{Float64}, X̃[:,lin_ind])
-        F = bunchkaufman(Hermitian(X̃' * Diagonal(w) * X̃), true)
+        F = bunchkaufman(Hermitian(X̃' * Diagonal(wts) * X̃), true)
     else
         lin_ind = collect(1:size(F, 2))
     end
-    β = F \ (X̃' * Diagonal(w) * y)
+    β = F \ (X̃' * Diagonal(wts) * y)
     Ψ = Hermitian(inv(F))
     ŷ = isempty(z) ? X̃ * β : hcat(X, z)[:,lin_ind] * β
-    X̃, y, β, Ψ, ŷ, w, lin_ind
+	X = X̃
+	wts = FrequencyWeights(collect(wts))
+	@pack! obj = X, y, β, Ψ, ŷ, wts
+	obj.vars = (obj.vars[1], obj.vars[2][lin_ind])
+	obj
 end
-"""
-    solve(estimator::ContinuousResponse,
-          X::AbstractMatrix{<:Number},
-          y::AbstractVector{<:Number},
-          z::AbstractVecOrMat{<:Number},
-          Z::AbstractMatrix{<:Number},
-          wts::AbstractVector)
-
-Solves continuous response models with potential features absorption.
-"""
-@views function solve(estimator::ContinuousResponse,
-                      X::AbstractMatrix{<:Number},
-                      y::AbstractVector{<:Number},
-                      z::AbstractVecOrMat{<:Number},
-                      Z::AbstractMatrix{<:Number},
-                      wts::AbstractVector)
+@views function fit!(obj::EconometricModel{<:ContinuousResponse})
+	@unpack estimator, X, y, z, Z, wts = obj
     if !isa(wts, Weights)
         wts = FrequencyWeights(wts)
     else
@@ -116,7 +150,11 @@ Solves continuous response models with potential features absorption.
     ŷ = isempty(z) ? X̃ * β : hcat(X, z)[:,lin_ind] * β
     û = y - ŷ
     ŷ .= y₀ .- û
-    X̃, y₀, β, Ψ, ŷ, w, lin_ind
+	X = X̃
+	y = y₀
+	@pack! obj = X, y, β, Ψ, ŷ
+	obj.vars = (obj.vars[1], obj.vars[2][lin_ind])
+	obj
 end
 """
     obtain_Ω(A::AbstractMatrix{<:Real},
@@ -142,25 +180,10 @@ Obtain Ω for a multinomial regression by building the matrix by blocks.
   end
   Hermitian(Σ)
 end
-"""
-    solve(estimator::NominalResponse,
-          X::AbstractMatrix{<:Number},
-          y::AbstractVector,
-          z::AbstractVecOrMat{<:Number},
-          Z::AbstractMatrix{<:Number},
-          wts::AbstractVector)
-
-Solves a multinomial logistic regression.
-"""
-@views function solve(estimator::NominalResponse,
-                      X::AbstractMatrix{<:Number},
-                      y::AbstractVector,
-                      z::AbstractVecOrMat{<:Number},
-                      Z::AbstractMatrix{<:Number},
-                      wts::AbstractVector)
+@views function fit!(obj::EconometricModel{<:NominalResponse})
+	@unpack estimator, X, y, z, Z, wts = obj
+	@unpack categories = estimator
     @assert isempty(z) && isempty(Z) "Nominal response models can only contain exogenous features"
-    @unpack categories = estimator
-    y = [ findfirst(isequal(x), categories) for x ∈ y ]
     b = mapreduce(elem -> (eachindex(categories) .== elem)', vcat, y)
     F = qr(X, Val(true))
     qrr = count(x -> abs(x) ≥ √eps(), diag(F.R))
@@ -209,7 +232,9 @@ Solves a multinomial logistic regression.
     Ψ = Hermitian(inv(bunchkaufman!(obtain_Ω(X, μ, wts))))
     ŷ = X * β
     β = collect(vec(β)[size(X, 2) + 1:end])
-    X, y, β, Ψ, ŷ, wts, 1:size(F, 2)
+	@pack! obj = X, y, β, Ψ, ŷ
+	obj.vars = (obj.vars[1], obj.vars[2][lin_ind])
+	obj
 end
 """
     solve(estimator::OrdinalResponse,
@@ -221,16 +246,10 @@ end
 
 Solves a proportional odds logistic regression.
 """
-@views function solve(estimator::OrdinalResponse,
-                      X::AbstractMatrix{<:Number},
-                      y::AbstractVector,
-                      z::AbstractVecOrMat{<:Number},
-                      Z::AbstractMatrix{<:Number},
-                      wts::AbstractVector)
+@views function fit!(obj::EconometricModel{<:OrdinalResponse})
+	@unpack estimator, X, y, z, Z, wts = obj
+	@unpack categories = estimator
     @assert isempty(z) && isempty(Z) "Ordinal response models can only contain exogenous features"
-    @unpack categories = estimator
-    y = [ findfirst(isequal(x), categories) for x ∈ y ]
-    @assert length(categories) > 2
     F = qr(X, Val(true))
     qrr = count(x -> abs(x) ≥ √eps(), diag(F.R))
     if qrr < size(F, 2)
@@ -292,5 +311,8 @@ Solves a proportional odds logistic regression.
     Ψ = Hermitian(A * Ψ * A')
     ŷ = X * β[1:size(X, 2)]
     β[ks:end] .= cumsum(vcat(β[ks], exp.(β[ks + 1:end])))
-    X, y, β, Ψ, ŷ, wts, collect(2:size(X, 2) + 1)
+	@pack! obj = X, y, β, Ψ, ŷ
+	# Fix the rank-deficient here
+	obj.vars = (obj.vars[1], obj.vars[2][collect(2:size(X, 2) + 1)])
+	obj
 end

src/statsbase.jl

@@ -42,7 +42,7 @@ function nullloglikelihood(obj::EconometricModel{<:Union{NominalResponse,Ordinal
 	sum(wᵢ * logpdf(Categorical(μ), yᵢ) for (yᵢ, wᵢ) ∈ zip(y, wts))
 end
 # StatsBase.score(obj::StatisticalModel) = error("score is not defined for $(typeof(obj)).")
-nobs(obj::EconometricModel) = sum(obj.w)
+nobs(obj::EconometricModel) = sum(weights(obj))
 dof(obj::EconometricModel{<:LinearModelEstimators}) =
 	length(coef(obj)) + dispersion(obj) + obj.iv
 dof(obj::EconometricModel{<:ContinuousResponse}) =
@@ -111,7 +111,7 @@ confint(obj::EconometricModel;
     σ * quantile(TDist(dof_residual(obj)), 1 - α / 2) |>
     (σ -> coef(obj) |>
         (β -> hcat(β .- σ, β .+ σ)))
-weights(obj::EconometricModel) = obj.w
+weights(obj::EconometricModel) = obj.wts
 isfitted(obj::EconometricModel) = !isempty(obj.Ψ)
 fitted(obj::EconometricModel) = obj.ŷ
 response(obj::EconometricModel) = obj.y

src/transformations.jl

@@ -74,7 +74,7 @@ Applies a transformation to a model component based on the model estimator.
 function transform end
 transform(estimator::LinearModelEstimators, wts::FrequencyWeights) = wts
 transform(estimator::BetweenEstimator, wts::FrequencyWeights) =
-    FrequencyWeights(Ones(length(estimator.groups)))
+    FrequencyWeights(ones(length(estimator.groups)))
 transform(estimator::ContinuousResponse,
           obj::AbstractVecOrMat{<:Number},
           wts::FrequencyWeights) = isempty(obj) ? obj : within(obj, estimator.groups, wts)