#=
Decomposition of Return
    rₜ=√hₜ zₜ, zₜ~N(0,1) so rₜ|Fₜ₋₁~N(0,hₜ). hₜ=Var(rₜ|Fₜ₋₁)

GARCH(p,q)(Bollerslev 1986)

    hₜ₊₁=ω+Σβⱼhₜ₋ⱼ+Σαᵢ(rₜ₋ᵢ²-hₜ₋ᵢ)

    αᵢ,βⱼ,ω≧0

EGARCH(p,q)(Nelson 1991)

    loghₜ₊₁=ω+Σαᵢg(zₜ)+Σβⱼloghₜ₋ⱼ

    where g(z)=θ₁z+θ₂(|z|-E|z|) (Leverage function)

REGARCH(1,1,K)(Hansen and Zhuo 2015)

    GARCH equation

        loghₜ₊₁=ω + τ(zₜ) + βloghₜ + γ′uₜ

        where τ(z)=τ₁z+τ₂(z²-1)

    Measurement Equation

        log xₖ,ₜ = ξₖ + ϕₖ log hₜ + δ₍ₖ₎(zₜ) + uₖ,ₜ for k=1,...,K (we assume ϕₖ=1)

        where δ₍ₖ₎(z)=δₖ,₁z + δₖ,₂(z²-1)

=#





module GARCHModels


    using Optim
    using Optim: TwiceDifferentiable, optimize, minimizer
    using FiniteDiff, ForwardDiff
    using LinearAlgebra, Statistics, Distributions
    using DataFrames

    export REGARCH, GARCH, EGARCH

    mutable struct Volatilityestimation
        parameter
        robustcov
        conditionalvariance
        LL
    end

    subscript(i::Integer) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)))

    function GARCH_coefnames(p,q)
        names = Array{String, 1}()
        push!(names,"ω")
        for k in 1:p
            push!(names, "β"*subscript(k))
        end
        for k in 1:q
            push!(names, "α"*subscript(k))
        end
        push!(names,"h0")
        return names
    end
    function GARCH_startingvals(p,q,rt)
        x0 = zeros(p+q+2)
        x0[1]=0.01
        x0[2:p+1] = 0.7/p
        x0[p+2:p+1+q] .= 0.2/q
        x0[end] = (var(rt))
        return x0
    end
    function GARCH_constraints(p,q)
        lower = -ones(p+q+2)
        upper = ones(p+q+2)
        lower[end] = 0
        return lower, upper
    end
    function GARCH_loglik(params::AbstractVector{T},p,q,rt,output) where T
    # Return data must be a Vector, not a Matrix.
        lower, upper = GARCH_constraints(p,q)
        all(lower.<params.<upper) || return -Inf
        n = size(rt)[1]
        presample = max(p,q)
        names = Array{String, 1}()

        ω = params[1]
        β = params[2:p+1]
        α = params[p+2:p+q+1]
        h0 = params[end]
        ht = zeros(T,n)
        zt = zeros(T,n)
        ht[1:presample] .= h0
        zt[1:presample] = rt[1:presample]./sqrt.(ht[1:presample])
        for t in presample+1:n
            mht = ω
            for i in 1:p
                mht += β[i]*ht[t-i]#log(xt[t,k]) - ξ[k] - log(ht[t]) - δ₁[k]*zt[t] - δ₂[k]*(zt[t]^2 - 1)
            end
            for j in 1:q
                mht += α[j]*(rt[t-j]^2 - ht[t-j])
            end
            ht[t] = mht
            zt[t] = rt[t]./sqrt.(ht[t])
        end
        LL = -log.(ht) .- zt.^2
        if output=="LL"
            return sum(LL)
        elseif output=="ht"
            return ht
        elseif output=="LLs"
            return LL
        end
    end
    function GARCH_scores(params,p,q,rt)
        f = x -> GARCH_loglik(x,p,q,rt,"LLs")
        S = ForwardDiff.jacobian(f, params)
    end
    function GARCH_InformationMatrix(params,p,q,rt)
        g = x -> GARCH_loglik(x,p,q,rt,"LL")
        H = ForwardDiff.hessian(g, params)
        J = -H
    end

    function GARCH_fit(p,q,rt,xt)
        obj = params -> -GARCH_loglik(params,p,q,rt,"LL")
        params = GARCH_startingvals(p,q,rt)
        res = optimize(obj, params, autodiff=:forward,iterations =10000,show_trace=true)
        params .= Optim.minimizer(res)
    end

    function GARCH(p,q,rt)
        params = GARCH_fit(p,q,rt)
        I = GARCH_InformationMatrix(params,p,q,rt)
        S = GARCH_scores(params,p,q,rt)
        ht = GARCH_loglik(params,p,q,rt,"ht")
        J = S'*S
        robustcov = I\J/I
        LL = GARCH_loglik(params,p,q,rt,"LL")
        coefnames = GARCH_coefnames(K)
        z = params./sqrt.(diag(robustcov))
        pvalue = (1 .- cdf.(Normal(),abs.(z)))*2
        parameters = DataFrame([coefnames,params,sqrt.(diag(robustcov)),z,pvalue])
        DataFrames.rename!(parameters,["Parameter","Estimate","Std","z","P(>|z|)"])
        res = Volatilityestimation(parameters,robustcov,ht,LL)
    end


    function EGARCH_coefnames(p,q)
        names = Array{String, 1}()
        push!(names,"ω")
        for k in 1:p
            push!(names, "β"*subscript(k))
        end
        for k in 1:q
            push!(names, "α"*subscript(k))
        end
        for i in 1:2
            push!(θ, "θ"*subscript(k))
        end
        push!(names,"logh0")
        return names
    end

    function EGARCH_startingvals(p,q,rt)
        x0 = zeros(p+q+4)
        x0[1]=0
        x0[2:p+1] = 0.7/p
        x0[p+2:p+1+q] .= 0.2/q
        x0[p+q+2:p+q+3] .= 0
        x0[end] = (var(rt))
        return x0
    end

    function EGARCH_constraints(p,q)
        lower = -ones(p+q+4)
        upper = ones(p+q+4)
        lower[end] = -10
        return lower, upper
    end

    function EGARCH_loglik(params::AbstractVector{T},p,q,rt,output) where T
    # Return data must be a Vector, not a Matrix.
        lower, upper = EGARCH_constraints(p,q)
        all(lower.<params.<upper) || return -Inf

        n = size(rt)[1]
        presample = max(p,q)
        names = Array{String, 1}()

        ω = params[1]
        β = params[2:p+1]
        α = params[p+2:p+q+1]
        θ = params[p+q+2:p+q+3]
        logh0 = params[end]
        ht = zeros(T,n)
        zt = zeros(T,n)
        ht[1:presample] .= exp(logh0)
        zt[1:presample] = rt[1:presample]./sqrt.(ht[1:presample])
        for t in presample+1:n
            mht = ω
            for i in 1:p
                mht += β[i]*log(ht[t-i])#log(xt[t,k]) - ξ[k] - log(ht[t]) - δ₁[k]*zt[t] - δ₂[k]*(zt[t]^2 - 1)
            end
            for j in 1:q
                mht += α[j]*(abs(zt[t-j]) - sqrt(2/pi))
            end
            ht[t] = mht
            zt[t] = rt[t]./sqrt.(ht[t])
        end
        LL = -log.(ht) .- zt.^2
        if output=="LL"
            return sum(LL)
        elseif output=="ht"
            return ht
        elseif output=="LLs"
            return LL
        end
    end
    function EGARCH_scores(params,p,q,rt)
        f = x -> EGARCH_loglik(x,p,q,rt,"LLs")
        S = ForwardDiff.jacobian(f, params)
    end
    function EGARCH_InformationMatrix(params,p,q,rt)
        g = x -> EGARCH_loglik(x,p,q,rt,"LL")
        H = ForwardDiff.hessian(g, params)
        J = -H
    end

    function EGARCH_fit(p,q,rt,xt)
        obj = params -> -EGARCH_loglik(params,p,q,rt,"LL")
        params = EGARCH_startingvals(p,q,rt)
        res = optimize(obj, params, autodiff=:forward,iterations =10000,show_trace=true)
        params .= Optim.minimizer(res)
    end

    function EGARCH(p,q,rt)
        params = EGARCH_fit(p,q,rt)
        I = EGARCH_InformationMatrix(params,p,q,rt)
        S = EGARCH_scores(params,p,q,rt)
        ht = EGARCH_loglik(params,p,q,rt,"ht")
        J = S'*S
        robustcov = I\J/I
        LL = EGARCH_loglik(params,p,q,rt,"LL")
        coefnames = EGARCH_coefnames(K)
        z = params./sqrt.(diag(robustcov))
        pvalue = (1 .- cdf.(Normal(),abs.(z)))*2
        parameters = DataFrame([coefnames,params,sqrt.(diag(robustcov)),z,pvalue])
        DataFrames.rename!(parameters,["Parameter","Estimate","Std","z","P(>|z|)"])
        res = Volatilityestimation(parameters,robustcov,ht,LL)
    end

    function REGARCH_coefnames(p,q)
        names = Array{String, 1}()
        push!(names,"ω")
        push!(names, "β")
        push!(names, "τ₁")
        push!(names, "τ₂")
        for k in 1:K
            push!(names, "γ"*subscript(k))
        end
        for k in 1:K
            push!(names, "ξ"*subscript(k))
        end
        for k in 1:K
            push!(names,"δ₁"*subscript(k))
        end
        for k in 1:K
            push!(names,"δ₂"*subscript(k))
        end
        push!(names,"logh0")
        return names
    end
    function REGARCH_startingvals(K,rt)
        x0 = zeros(5+4*K)
        x0[1]=0.01
        x0[2] = 0.8
        x0[5:4+K] .= 0.3/K
        x0[5+K:4+2*K] .= -0.3
        x0[end] = log(var(rt))
        return x0
    end
    function REGARCH_constraints(K)
        lower = -ones(4*K+5)
        upper = ones(4*K+5)
        upper[5+2*K:4+3*K] .= 2
        lower[end] = -50
        return lower, upper
    end
    function REGARCH_loglik(params::AbstractVector{T},K,rt,xt,output) where T
    # Return data must be a Vector, not a Matrix.
        size(rt)[1]==size(xt)[1] || return error("Check the length of the data")
        lower, upper = REGARCH_constraints(K)
        all(lower.<params.<upper) || return -Inf
        n = size(rt)[1]
        presample = 1
        nparam = 5+5*K
        ω = params[1]
        β = params[2]
        τ₁ = params[3]
        τ₂ = params[4]
        γ = params[5:4+K]
        ξ = params[5+K:4+2*K]
        δ₁ = params[5+2*K:4+3*K]
        δ₂ = params[5+3*K:4+4*K]
        logh0 = params[end]
        ht = zeros(T,n)
        ut = zeros(T,n,K)
        zt = zeros(T,n)
        ht[1:presample] .= exp(logh0)
        ut[1:presample] .= 0
        zt[1:presample] = rt[1:presample]./sqrt.(ht[1:presample])
        for t in presample+1:n
            mlht = ω + β*log(ht[t-1]) .+ τ₁*zt[t-1] .+ τ₂*(zt[t-1]^2-1) .+ sum(γ.*(ut[t-1,:]))
            ht[t] = exp(mlht)
            zt[t] = rt[t]./sqrt.(ht[t])
            for k in 1:K
                ut[t,k] = log(xt[t,k]) - ξ[k] - log(ht[t]) - δ₁[k]*zt[t] - δ₂[k]*(zt[t]^2 - 1)
            end
        end
        Σ = ut'*ut/n
        LL = -log.(ht) .- zt.^2 .- log(det(Σ))
        if output=="LL"
            return sum(LL)
        elseif output=="ht"
            return ht
        elseif output=="LLs"
            return LL
        end
    end
    function REGARCH_scores(params,K,rt,xt)
        f = x -> REGARCH_loglik(x,K,rt,xt,"LLs")
        S = ForwardDiff.jacobian(f, params)
    end
    function REGARCH_InformationMatrix(params,K,rt,xt)
        g = x -> REGARCH_loglik(x,K,rt,xt,"LL")
        H = ForwardDiff.hessian(g, params)
        J = -H
    end

    function REGARCH_fit(K,rt,xt)
        obj = params -> -REGARCH_loglik(params,K,rt,xt,"LL")
        params = REGARCH_startingvals(K,rt)
        res = optimize(obj, params, autodiff=:forward,iterations =10000,show_trace=true)
        params .= Optim.minimizer(res)
    end

    function REGARCH(K,rt,xt)
        params = REGARCH_fit(K,rt,xt)
        I = REGARCH_InformationMatrix(params,K,rt,xt)
        S = REGARCH_scores(params,K,rt,xt)
        ht = REGARCH_loglik(params,K,rt,xt,"ht")
        J = S'*S
        robustcov = I\J/I
        LL = REGARCH_loglik(params,K,rt,xt,"LL")
        coefnames = REGARCH_coefnames(K)
        z = params./sqrt.(diag(robustcov))
        pvalue = (1 .- cdf.(Normal(),abs.(z)))*2
        parameters = DataFrame([coefnames,params,sqrt.(diag(robustcov)),z,pvalue])
        DataFrames.rename!(parameters,["Parameter","Estimate","Std","z","P(>|z|)"])
        res = Volatilityestimation(parameters,robustcov,ht,LL)
    end



end
