GARCH in Turing.jl

garch.jl

using Turing
using CSVFiles
using DataFrames
using Dates
using StatsPlots

# The TV syntax allows sampling to be type-stable.
@model garch(r, ::Type{TV}=Vector{Float64}) where {TV} = begin
    T = length(r)

    # Declare priors on parameters.
    omega ~ TruncatedNormal(0,1,0,1)
    alpha ~ TruncatedNormal(0,1,0,1)
    beta ~ TruncatedNormal(0,1,0,1)

    # Initialize conditional variance vector.
    h = TV(undef, T)

    # Treat the first h as a parameter to estimate.
    h[1] ~ TruncatedNormal(0,1,0,1)

    # Observe each data point.
    for t in 2:T
        h[t] = omega + alpha * r[t-1]^2 + beta * h[t-1]
		r[t] ~ Normal(0, sqrt(h[t]))
    end
end

# Load S&P 500 prices.
df = DataFrame(load(joinpath(@__DIR__, "snp-prices.csv")))

# Add dates.
date_format = dateformat"Y-m-d"
date_str = string.(df.DATE)
df.DATE = [Date(string(d[1:4], "-", d[5:6], "-", d[7:8]), date_format) for d in date_str]

# Make a returns vector.
returns = [(df.spindx[i] - df.spindx[i-1])/df.spindx[i-1] for i in 2:size(df, 1)]

# Use the past 600 daily returns.
r = returns[(end-600):end]

# Demean the series.
r = r .- mean(r)

chain = sample(garch(r), NUTS(5_000, 0.65))

function drawpath(alpha, beta, omega, h1, r)
    h = Vector(undef, length(r))
    h[1] = h1
    T = length(r)

    for t in 2:T
        h[t] = omega + alpha*r[t-1]^2 + beta*h[t-1]
    end

    return h
end

function simulate_h(chain, r)
    theta = get_params(chain)
    draws = Vector{Vector{Float64}}(undef, length(chain))
    T = length(chain)

    for i in 1:length(chain)
        draws[i] = drawpath(theta.alpha[i], theta.beta[i], theta.omega[i], theta.h[i], r)
    end

    return draws
end

# Save the chain plot.
chain_plot = plot(chain, size=(800,600))
savefig(chain_plot, "chainplot.pdf")

df_one = DataFrame(Date=df.DATE[(end-600):end], Returns=r)
# p1 = @df df_one plot(:Date, :Volatility, title="Conditional Volatility")
# p2 = @df df_one plot(:Date, :Returns, title="Returns")
# p = plot(p1, p2, layout=(2,1), legend=false)
# savefig(p, "plot1.pdf")

draws = simulate_h(chain, r)
# draws = draw_h(mean(chain, :alpha), mean(chain, :beta), mean(chain, :omega), mean(chain, :h), length(r))
draws_plot_paths = plot(df_one.Date, rand(draws, 100), 
    legend=false, 
    alpha=0.3, 
    color_palette=:blues, 
    title="Conditional Volatility Paths",
    dpi=350);
draws_plot_returns = plot(df_one.Date, df_one.Returns, 
    legend=false, 
    title="Returns",
    xticks=false, 
    dpi=350);
draws_plot = plot(draws_plot_returns, draws_plot_paths, layout=(2,1), xrotation=15, dpi=350);
savefig(draws_plot, "drawsplot-all.png")

mu = mean(draws)
sd = std(draws)

sds = [i .* sd for i in 1:1:3]


p = plot();
al = 0.60
for sd in sds
    plot!(p, df_one.Date, mu, ribbon=sd, fillalpha=al, legend=false,
        linecolor=:transparent, color=:steelblue, dpi=350, title="Conditional Volatility Posterior");
    global al = al / 2
end
plot!(p, df_one.Date, mu, dpi=350, color=:blue);

savefig(plot(draws_plot_returns, p, layout=(2,1)), "drawpaths-region.png")
# Write the tables
open("summary_table1.txt", "w") do f
    m = summarystats(chain).df
    display(TextDisplay(f), "text/latex", m)
end

open("summary_table2.txt", "w") do f
    m = quantile(chain).df
    display(TextDisplay(f), "text/latex", m)
end

# The TV syntax allows sampling to be type-stable.