Upgrade to the GARCH implementation

REQUIRE

@@ -3,3 +3,4 @@ Distributions
 StatsBase
 TimeSeries 0.5.11
 Optim
+NLopt

src/GARCH.jl

@@ -2,120 +2,115 @@
 # Copyright 2013 Andrey Kolev
 # Distributed under MIT license (see LICENSE.md)
 
-type GarchFit
-  data::Vector
-  params::Vector
-  llh::Float64
-  status::Symbol
-  converged::Bool
-  sigma::Vector
-  hessian::Array{Float64,2}
-  cvar::Array{Float64,2}
-  secoef::Vector
-  tval::Vector
+module GARCH
+using NLopt, Distributions
+
+export fit_GARCH, predict
+
+immutable GarchFit
+    ɛ::Vector{Float64}
+    ɛ²::Vector{Float64}
+    params::Vector{Float64}
+    llh::Float64
+    status::Symbol
+    converged::Bool
+    sigma::Vector{Float64}
+    hessian::Array{Float64,2}
+    cvar::Array{Float64,2}
+    secoef::Vector{Float64}
+    tval::Vector{Float64}
+end
+
+function Base.show(io::IO, fit::GarchFit)
+    pnorm(x) = 0.5 * (1 + erf(x / sqrt(2)))
+    prt(x) = 2 * (1 - pnorm(abs(x)))
+    jbstat, jbp = jbtest(fit.ɛ ./ fit.sigma)
+
+    @printf io "Fitted garch model \n"
+    @printf io " * Coefficient(s):    %-15s%-15s%-15s\n" "α₀" "α₁" "β₁"
+    @printf io "%-22s%-15.5g%-15.5g%-15.5g\n" "" fit.params[1] fit.params[2] fit.params[3]
+    @printf io " * Log Likelihood: %.5g\n" fit.llh
+    @printf io " * Converged: %s\n" fit.converged
+    @printf io " * Solver status: %s\n\n" fit.status
+    @printf io " * Standardised Residuals Tests:\n"
+    @printf io "   %-26s%-15s%-15s\n" "" "Statistic" "tp-Value"
+    @printf io "   %-21s%-5s%-15.5g%-15.5g\n\n" "Jarque-Bera Test" "χ²" jbstat jbp
+    @printf io " * Error Analysis:\n"
+    @printf io "   %-7s%-15s%-15s%-15s%-15s\n" "" "Estimate" "Std.Error" "t value" "Pr(>|t|)"
+    @printf io "   %-7s%-15.5g%-15.5g%-15.5g%-15.5g\n" "α₀" fit.params[1] fit.secoef[1] fit.tval[1] prt(fit.tval[1])
+    @printf io "   %-7s%-15.5g%-15.5g%-15.5g%-15.5g\n" "α₁" fit.params[2] fit.secoef[2] fit.tval[2] prt(fit.tval[2])
+    @printf io "   %-7s%-15.5g%-15.5g%-15.5g%-15.5g\n" "β₁" fit.params[3] fit.secoef[3] fit.tval[3] prt(fit.tval[3])
 end
 
-function Base.show(io::IO ,fit::GarchFit)
-  pnorm(x) = 0.5*(1+erf(x/sqrt(2)))
-  prt(x) = 2*(1-pnorm(abs(x)))
-  @printf io "Fitted garch model \n"
-  @printf io " * Coefficient(s): \tomega \t\talpha \t\tbeta\n"
-  @printf io "   \t\t\t%f\t%f\t%f\n" fit.params[1] fit.params[2] fit.params[3]
-  @printf io " * Log Likelihood: %f\n" fit.llh
-  @printf io " * Converged: %s\n" fit.converged
-  @printf io " * Solver status: %s\n\n" fit.status
-  println(io," * Standardised Residuals Tests:")
-  println(io,"   \t\t\t\tStatistic\tp-Value")
-  jbstat,jbp = jbtest(fit.data./fit.sigma);
-  @printf io "   Jarque-Bera Test\t\U1D6D8\u00B2\t%.6f\t%.6f\n\n" jbstat jbp
-  println(io," * Error Analysis:")
-  println(io,"   \t\tEstimate\t\Std.Error\tt value \tPr(>|t|)")
-  @printf io "   omega\t%f\t%f\t%f\t%f\n" fit.params[1] fit.secoef[1] fit.tval[1] prt(fit.tval[1])
-  @printf io "   alpha\t%f\t%f\t%f\t%f\n" fit.params[2] fit.secoef[2] fit.tval[2] prt(fit.tval[2])
-  @printf io "   beta \t%f\t%f\t%f\t%f\n"  fit.params[3] fit.secoef[3] fit.tval[3] prt(fit.tval[3])
+function cdHessian(par, LLH)
+    eps = 1e-4 * par
+    n = length(par)
+    H = zeros(n, n)
+    for i in 1:n
+        for j in 1:n
+            x₁ = copy(par)
+            x₁[i] += eps[i]
+            x₁[j] += eps[j]
+            x₂ = copy(par)
+            x₂[i] += eps[i]
+            x₂[j] -= eps[j]
+            x₃ = copy(par)
+            x₃[i] -= eps[i]
+            x₃[j] += eps[j]
+            x₄ = copy(par)
+            x₄[i] -= eps[i]
+            x₄[j] -= eps[j]
+            H[i,j] = (LLH(x₁) - LLH(x₂) - LLH(x₃) + LLH(x₄)) / (4 .* eps[i] * eps[j])
+        end
+    end
+    return H
 end
 
-function cdHessian(par,LLH)
-  eps = 1e-4 * par
-  n = length(par)
-  H = zeros(n,n)
-  for(i = 1:n)
-    for(j = 1:n)
-      x1 = copy(par) 
-      x1[i] += eps[i]
-      x1[j] += eps[j] 
-      x2 = copy(par)
-      x2[i] += eps[i]
-      x2[j] -= eps[j]
-      x3 = copy(par)
-      x3[i] -= eps[i]
-      x3[j] += eps[j]
-      x4 = copy(par)
-      x4[i] -= eps[i]
-      x4[j] -= eps[j]
-      H[i,j] = (LLH(x1)-LLH(x2)-LLH(x3)+LLH(x4)) / (4.*eps[i]*eps[j])
+
+function calculate_volatility_process(ɛ²::Vector, α₀, α₁, β₁)
+    h = similar(ɛ²)
+    h[1] = mean(ɛ²)
+    for i = 2:length(ɛ²)
+        h[i] = α₀ + α₁ * ɛ²[i - 1] + β₁ * h[i - 1]
     end
-  end
-  H
+    return h
 end
 
-function garchLLH(rets::Vector,x::Vector)
-  rets2   = rets.^2;
-  T = length(rets); 
-  ht = zeros(T);
-  omega,alpha,beta = x;
-  ht[1] = sum(rets2)/T;
-  for i=2:T
-    ht[i] = omega + alpha*rets2[i-1] + beta * ht[i-1];
-  end
-  -0.5*(T-1)*log(2*pi)-0.5*sum( log(ht) + (rets./sqrt(ht)).^2 );
+
+function fit_GARCH_LLH(ɛ²::Vector, x::Vector)
+    T = length(ɛ²)
+    α₀, α₁, β₁ = x
+    h = calculate_volatility_process(ɛ², α₀, α₁, β₁)
+    return -0.5 * (T - 1) * log(2π) - 0.5 * sum(log(h) .+ ɛ² ./ h)
 end
 
 function predict(fit::GarchFit)
- omega, alpha, beta = fit.params;
- rets = fit.data
- rets2   = rets.^2;
- T = length(rets); 
- ht    = zeros(T);
- ht[1] = sum(rets2)/T;
- for i=2:T
-    ht[i] = omega + alpha*rets2[i-1] + beta * ht[i-1];
- end
- sqrt(omega + alpha*rets2[end] + beta*ht[end]);
+    α₀, α₁, β₁ = fit.params
+    ɛ² = fit.ɛ²
+    h = calculate_volatility_process(ɛ², α₀, α₁, β₁)
+    return sqrt(α₀ + α₁ * ɛ²[end] + β₁ * h[end])
 end
 
-function garchFit(data::Vector)
-  rets = data
-  rets2   = rets.^2;
-  T = length(rets); 
-  ht = zeros(T);
-  function garchLike(x::Vector, grad::Vector)
-    omega,alpha,beta = x;
-    ht[1] = sum(rets2)/T;
-    for i=2:T
-      ht[i] = omega + alpha*rets2[i-1] + beta * ht[i-1];
+function fit_GARCH(ɛ::Vector)
+    ɛ² = ɛ .^ 2
+    T = length(ɛ²)
+    h = zeros(T)
+    function fit_GARCH_like(x::Vector, grad::Vector)
+        α₀, α₁, β₁ = x
+        h = calculate_volatility_process(ɛ², α₀, α₁, β₁)
+        sum(log(h) .+ ɛ² ./ h)
     end
-    sum( log(ht) + (rets./sqrt(ht)).^2 );
-  end
-  opt = Opt(:LN_SBPLX,3)
-  lower_bounds!(opt,[1e-10, 0.0, 0.0])
-  upper_bounds!(opt,[1; 0.3; 0.99])
-  min_objective!(opt, garchLike)
-  (minf,minx,ret) = Optim.optimize(opt, [1e-5, 0.09, 0.89])
-  converged = minx[1]>0 && all(minx[2:3].>=0) && sum(minx[2:3])<1.0
-  H = cdHessian(minx,x->garchLLH(rets,x))
-  cvar = -inv(H)
-  secoef = sqrt(diag(cvar))
-  tval = minx./secoef
-  out = GarchFit(data, minx, -0.5*(T-1)*log(2*pi)-0.5*minf, ret, converged, sqrt(ht),H,cvar,secoef,tval)
+    opt = Opt(:LN_SBPLX, 3)
+    lower_bounds!(opt, [1e-10, 0.0, 0.0])
+    upper_bounds!(opt, [1, 0.3, 0.99])
+    min_objective!(opt, fit_GARCH_like)
+    (minf, minx, ret) = NLopt.optimize(opt, [1e-5, 0.09, 0.89])
+    converged = minx[1] > 0 && all(minx[2:3] .>= 0) && sum(minx[2:3]) < 1.0
+    H = cdHessian(minx, x -> fit_GARCH_LLH(ɛ², x))
+    cvar = -inv(H)
+    secoef = sqrt(diag(cvar))
+    tval = minx ./ secoef
+    return GarchFit(ɛ, ɛ², minx, -0.5 * (T - 1) * log(2π) - 0.5 * minf, ret, converged, sqrt(h), H, cvar, secoef, tval)
 end
 
-# function garchPkgTest()
-#   println("Running GARCH package test...")
-#   try
-#     include(Pkg.dir("GARCH", "test","GARCHtest.jl"))
-#     println("All tests passed!")
-#   catch err
-#     throw(err)
-#   end
-# end
+end #GARCH

src/TimeModels.jl

@@ -3,26 +3,23 @@ module TimeModels
 using Base.Dates
 using Distributions
 using StatsBase
-using TimeSeries 
+using TimeSeries
 using Optim
 
 import Base: show
 
 export
   # Kalman exports
   StateSpaceModel,
-  KalmanFiltered, 
+  KalmanFiltered,
   KalmanSmoothed,
   simulate,
   kalman_filter,
   kalman_smooth,
-  fit, 
+  fit,
   # ARIMA exports
   arima_statespace,
   arima,
-  # GARCH exports
-  garchFit,
-  predict,
   # diagnostic tests exports
   jbtest
 
@@ -39,4 +36,4 @@ include("GARCH.jl")
 # Tests
 include("diagnostic_tests.jl")
 
-end 
+end

test/GARCH.jl

@@ -0,0 +1,14 @@
+using FactCheck
+using TimeModels
+using TimeModels.GARCH
+
+facts("GARCH") do
+    context("Consistency with R's fGarch") do
+        filename = Pkg.dir("TimeModels", "test", "data", "random process.csv")
+        ret = readcsv(filename)[:, 1]
+        ret = ret .- mean(ret) #Note: GARCH does not remove the mean, so we have to remove it ourselves.
+        fit = fit_GARCH(ret)
+        param = [2.469347e-06, 1.142268e-01, 8.691734e-01] #R fGarch garch(1,1) estimated params on the zeroized return.
+        @fact fit.params --> roughly(param, atol=1e-3)
+    end
+end

test/all.jl

@@ -1,6 +1,10 @@
 # this file is useful for interactive REPL testing
 
+module TestCapsule
+
 # include("arima.jl")
 # include("diagnostic_tests.jl")
-# include("garch.jl") 
+include("garch.jl")
 include("kalman_filter.jl")
+
+end #TestCapsule