Recursive GARCH

garchrec.jl

function ht(theta, Y, t)
    if t == 1
        return theta[2]^2 / (1 - theta[3]^2 - theta[4]^2)
    end

    theta[2]^2 + theta[3]^2 * (Y[t-1] - theta[1])^2 + theta[4]^2 * ht(theta, Y, t-1)
end
    
function score(theta, Y, t)
    var = ht(theta, Y, t)
    if var <= 0
        return oftype(var, Inf)
    end
    d = Normal(theta[1], sqrt(var))
    logpdf(d, Y[t])
end

function scores(theta, Y)
    T = length(Y)
    map(x -> score(theta, Y, x), 1:T)
end
    
function loglik(theta, y)
    -sum(scores(theta, y))
end

@time optimum = Optim.optimize(x -> loglik(x, R), X0, LBFGS())
XHat = optimum.minimum

#225.945721 seconds (9.69 M allocations: 227.654 MB, 0.02% gc time)
#Out[13]:
#4-element Array{Float64,1}:
# 0.0464168
# 0.119154 
# 0.276747 
# 0.954312