Add negative binomial distribution (#99)

* added fitted quantiles and SDMforecast

* update default verbose

* remove fitted quantiles

* update

* add ane to test folder

* add NegativeBinomial distribution

* add negbin to docs

* fix special functions dep

Author: guilhermebodin

Project.toml

@@ -1,7 +1,7 @@
 name = "ScoreDrivenModels"
 uuid = "4a87933e-d659-11e9-0e65-7f40dedd4a3a"
 authors = ["guilhermebodin <[email protected]>, raphaelsaavedra <[email protected]>"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -12,7 +12,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 julia = "1"
-Distributions = "0.21"
+Distributions = "0.23"
 Optim = "0.20"
 SpecialFunctions = "0.8"
 

docs/src/manual.md

@@ -99,6 +99,7 @@ this is handled internally.
 |[`Gamma`](@ref)| ✓ | ✓  |
 |[`LogitNormal`](@ref)| ✓ | ✓  |
 |[`LogNormal`](@ref)| ✓ | ✓  |
+|[`NegativeBinomial`](@ref)| ✓ |  x  |
 |[`Normal`](@ref)| ✓ | ✓  |
 |[`Poisson`](@ref)| ✓ | ✓  |
 |[`TDist`](@ref)| ✓ | ✓  |
@@ -115,6 +116,7 @@ ScoreDrivenModels.Exponential
 ScoreDrivenModels.Gamma
 ScoreDrivenModels.LogitNormal
 ScoreDrivenModels.LogNormal
+ScoreDrivenModels.NegativeBinomial
 ScoreDrivenModels.Normal
 ScoreDrivenModels.Poisson
 ScoreDrivenModels.TDist

examples/ane.jl

@@ -0,0 +1,26 @@
+using ScoreDrivenModels, Plots
+
+# Convert data to vector
+y = Vector{Float64}(vec(readdlm("./test/data/ena_northeastern.csv")'))
+y_train = y[1:400]
+y_test = y[401:end]
+
+# Specify model: here we use lag 1 for trend characterization and lag 12 for seasonality characterization
+gas = Model([1, 2, 11, 12], [1, 2, 11, 12], LogNormal, 0.0; time_varying_params = [1])
+
+# Define initial_params with
+initial_params = dynamic_initial_params(y_train, gas)
+
+# Estimate the model via MLE
+f = fit!(gas, y_train; initial_params = initial_params)
+
+# Compare observations and in-sample estimates
+plot(y_train, label = "In-sample ANE")
+
+# Forecasts with 95% confidence interval
+forecast = forecast_quantiles(y_train, gas, 60; S = 1_000, initial_params = initial_params)
+
+plot(forecast.scenarios, color = "grey", width = 0.05, label = "")
+plot!(y[360:460], label = "ANE", color = "black", xlabel = "Months", ylabel = "GWmed", legend = :topright)
+plot!(forecast.quantiles, label = ["Quantiles" "" ""], color = "red", line = :dash)
+

examples/inflow_lognormal.jl

@@ -39,6 +39,7 @@ include("distributions/exponential.jl")
 include("distributions/gamma.jl")
 include("distributions/logitnormal.jl")
 include("distributions/lognormal.jl")
+include("distributions/negativebinomial.jl")
 include("distributions/normal.jl")
 include("distributions/poisson.jl")
 include("distributions/tdist.jl")

src/ScoreDrivenModels.jl

@@ -27,7 +27,7 @@ end
 function log_likelihood(::Type{Chi}, y::Vector{T}, param::Matrix{T}, n::Int) where T
     loglik = 0.0
     for t in 1:n
-        loglik += + (1 - param[t, 1]/2)*log(2) + (param[t, 1] - 1) * log(y[t]) - (y[t]^2)/2 - loggamma(param[t, 1]/2)
+        loglik += (1 - param[t, 1] / 2) * log(2) + (param[t, 1] - 1) * log(y[t]) - (y[t]^2)/2 - loggamma(param[t, 1]/2)
     end
     return -loglik
 end

src/distributions/chi.jl

@@ -8,6 +8,7 @@ const DISTS = [
     Gamma,
     LogitNormal,
     LogNormal,
+    NegativeBinomial,
     Normal,
     Poisson,
     TDist,
@@ -23,6 +24,7 @@ export Beta,
     Gamma,
     LogitNormal,
     LogNormal,
+    NegativeBinomial,
     Normal,
     Poisson,
     TDist,

src/distributions/common_interface.jl

@@ -0,0 +1,64 @@
+"""
+    NegativeBinomial
+
+* Parametrization
+parametrized in r, p
+
+* Score
+
+* Fisher Information
+
+* `time_varying_params` map.
+
+* Default link
+"""
+NegativeBinomial
+
+function score!(score_til::Matrix{T}, y::Int, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
+    score_til[t, 1] = digamma(y + param[t, 1]) - digamma(param[t, 1]) + log(param[t, 2])
+    score_til[t, 2] = param[t, 1] / param[t, 2] - y / (one(T) - param[t, 2])
+    return
+end
+
+function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
+    return error("Fisher information not implemented for NegativeBinomial distribution.")
+end
+
+function log_likelihood(::Type{NegativeBinomial}, y::Vector{Int}, param::Matrix{T}, n::Int) where T
+    loglik = zero(T)
+    for t in 1:n
+        loglik += loggamma(y[t] + param[t, 1]) - logfactorial(y[t]) - loggamma(param[t, 1]) + 
+                param[t, 1] * log(param[t, 2]) + y[t] * log(1 - param[t, 2])
+    end
+    return -loglik
+end
+
+# Links
+function link!(param_tilde::Matrix{T}, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T 
+    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
+    param_tilde[t, 2] = link(LogitLink, param[t, 2], zero(T), one(T))
+    return
+end
+function unlink!(param::Matrix{T}, ::Type{NegativeBinomial}, param_tilde::Matrix{T}, t::Int) where T 
+    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
+    param[t, 2] = unlink(LogitLink, param_tilde[t, 2], zero(T), one(T))
+    return
+end
+function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T 
+    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
+    aux.jac[2] = jacobian_link(LogitLink, param[t, 2], zero(T), one(T))
+    return
+end
+
+# utils
+function update_dist(::Type{NegativeBinomial}, param::Matrix{T}, t::Int) where T
+    return NegativeBinomial(param[t, 1], param[t, 2])
+end 
+
+function params_sdm(d::NegativeBinomial)
+    return Distributions.params(d)
+end
+
+function num_params(::Type{NegativeBinomial})
+    return 2
+end