Skip to content

Instantly share code, notes, and snippets.

@cpfiffer
Created November 26, 2019 14:32
Show Gist options
  • Save cpfiffer/92f3f1ebf886ac58f4365f6d07d214ad to your computer and use it in GitHub Desktop.
Save cpfiffer/92f3f1ebf886ac58f4365f6d07d214ad to your computer and use it in GitHub Desktop.
GARCH in Turing.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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment