garch.jl

function score(theta::Array, y::Array)
    T = size(y, 1)
    h = zeros(T)
    eps = y .- theta[1]
    h[1] = theta[2]^2 / (1 - theta[3]^2 - theta[4]^2)
    if h[1] < 0
        return -oftype(h[1], Inf)
    end
    for  t=2:T
        h[t] = theta[2]^2 + theta[3]^2 * eps[t-1]^2 + theta[4]^2 * h[t-1]
    end        
    dists = map(x -> Normal(theta[1], sqrt(x)), h)

    map((x, y) -> logpdf(x, y), dists, y)
end

function loglik(theta::Array, y::Array)
    -sum(score(theta, y))
end

optimum = Optim.optimize(x -> loglik(x, R), X0, LBFGS())
XHat = optimum.minimum
# 4-element Array{Float64,1}:
#  0.0464168
#  0.119154 
#  0.276747 
#  0.954312