reproducing Heston example for microstructure noise paper

example_heston.jl

using MicrostructureNoise, Distributions 
using Random, Dates, LinearAlgebra, DelimitedFiles
using Bridge
using StaticArrays
using DataFrames
using CSV

const R2 = SVector{2,Float64};

Random.seed!(12)

NBINS = 100

############################## Generate data from forward model ##############################
forwardmodel= [:Heston, :IGlimit][1]  # in case IGlimit, then we generate the volatility according to the derived limit behaviour of the IGMC

n = 4_000       # nr of observations
ρ = 0.0#-0.6        # corr between Wiener processes
η0 = 1e-6       # variance of extrinsic noise

# set pars in Heston model
μ = 0.05
κ = 7.0
θ = 0.04
σ = 0.6

#---------------------------------------------

struct Heston <: ContinuousTimeProcess{R2}
    μ::Float64
    κ::Float64
    θ::Float64
    σ::Float64
end

Bridge.b(s, x, P::Heston) = R2(P.μ*x[1], P.κ*(P.θ - x[2]))
Bridge.σ(s, x, P::Heston) = Diagonal(R2(sqrt(x[2])*x[1], P.σ*sqrt(x[2])))

#---------------------------------------------
struct IGlimit <: ContinuousTimeProcess{R2} 
    μ::Float64
    σ::Float64
end

Bridge.b(s, x, P::IGlimit) = R2(P.μ *x[1], 0.0)
Bridge.σ(s, x, P::IGlimit) = Diagonal(R2(sqrt(x[2])*x[1], sqrt(2)*x[2]))

#---------------------------------------------
T = 1.0

if forwardmodel== :Heston
    P = Heston(μ, κ, θ, σ)
else
    P = IGlimit(μ, σ)
end

nf = 10 # generate the process on a finer grid than it is observed

tfine0 = T .* sort(rand(n*nf-1))
tfine = [0.0; tfine0; T]

is = [i <= n-1 for i in 1:length(tfine)-2]
is = [true; is[randperm(length(is))]; true]
t = tfine[is];

u = R2(1., θ)  # initial value of diffusion process
W = sample(tfine, Wiener{R2}())

map!(v->R2(v[1], ρ*v[1] + sqrt(1-ρ^2)v[2]), W.yy, W.yy) # correlate Brownian motions

#----------------------------------------------
Xfine = solve(EulerMaruyama(), u, W, P)

xtrue = log.(first.(Xfine.yy[is]))
s0(t) = sqrt(Xfine.yy[searchsortedfirst(Xfine.tt, t)][2])
y = copy(xtrue) 

y .+= randn(n + 1) * sqrt(η0);  # observe X_t = log S_t with extrinsic noise


############################## Perform inference using MicrostructureNoise. ##############################


α = 5.0 # Initial smoothness hyperparameter guess
σα = 3.0 # Random walk step size for smoothness hyperparameter


prior = MicrostructureNoise.Prior(
N = NBINS,

α1 = 0.10,
β1 = 0.10,

αη = 0.01, 
βη = 0.01,

Πα = LogNormal(1., 2.5),
μ0 = 0.0,
C0 = 5.0
);

td, θs, ηs, αs, p = MicrostructureNoise.MCMC(prior, t, y, α, σα, 3_000, printiter=500);
# θs: iterates of θ
# ηs: iterates of η
# αs: iterates of α

println("acceptance probability $p")

# further processing of iterates
posterior = MicrostructureNoise.posterior_volatility(td, θs)


# Julia plotting
tt, mm = MicrostructureNoise.piecewise(posterior.post_t, posterior.post_median[:])

using Plots

Plots.plot(tt, mm, label="posterior median", linewidth=0.5, color="dark blue")
Plots.plot!(t[1:10:end], (s0.(t[1:10:end])).^2, label="true volatility", color="red")


plot!(MicrostructureNoise.piecewise(posterior.post_t, posterior.post_qlow[:])...,
    fillrange = MicrostructureNoise.piecewise(posterior.post_t, posterior.post_qup[:])[2], 
    fillalpha = 0.2,
    alpha = 0.0,
    fillcolor="darkblue",
    title="Volatility estimate", label="marginal $(round(Int,posterior.qu*100))% credibility band")

# display(p)

# exporting data to csv 

# posterior has fieldnames
#(:post_t, :post_qlow, :post_median, :post_qup, :post_mean, :post_mean_root, :qu)

# write to csv files
pp = posterior 

d = DataFrame(bin = 1:length(pp.post_qlow), t=pp.post_t[2:end], post_qlow= pp.post_qlow, post_qup = pp.post_qup, post_mean= pp.post_mean, post_mean_root=pp.post_mean_root, post_median = pp.post_median)
s0_df = DataFrame(t = t, s0 =s0.(t) )

naming = string(forwardmodel)*"_"*string(NBINS)*".csv"
CSV.write("/Users/frankvandermeulen/Sync/DOCUMENTS/onderzoek/code/microstructure/out/posteriorinfo"*naming,  d)
CSV.write("/Users/frankvandermeulen/Sync/DOCUMENTS/onderzoek/code/microstructure/out/truevol"*naming,  s0_df)


# R ggplot plotting
using RCall

@rput d
@rput s0_df

R"""
library(tidyverse)
L <- nrow(d)
breaks=c(0,d$t)
d <-  d %>%  mutate(time=bin,xmin = breaks[-(L+1)],xmax=breaks[-1])%>%
  mutate(post_qlow=sqrt(post_qlow),post_qup=sqrt(post_qup))

  d_step <- data.frame(x=breaks,y=c(d$post_mean_root,d$post_mean_root[L]))

  ggplot() +
  geom_rect(data=d,aes(xmin=xmin,xmax=xmax,ymin = post_qlow, ymax = post_qup), fill = "lightsteelblue1")+
  geom_step(data=d_step, aes(x=x,y=y),colour='black',size=1,alpha=1)+ 
  geom_line(data=s0_df, aes(x=t, y=s0), colour='red',size=.5,alpha=0.8)+
  ylab("")+xlab("") #+ggtitle('Posterior mean  and 95% marginal credible bands')

"""