In [ ]:
### The ou_bm.ipynb IJulia notebook tests the OU-BM inferential tools in RDE
### The tests are based on numerical and visual output
### They overlap with but also complement the tests in ou_bm.jl
### Load the packages needed for running notebook simulations
using Brownian, Gadfly, RDE
In [ ]:
### Various configurations
#Configure the plotting environment of the Gadfly package
set_default_plot_size(18cm, 8cm);
In [ ]:
### Simulation 1: comparing exact and approx MLE estimators of OU-BM parameters
# Specify Brownian motion
t = 0:1/2^2:5;
bm = BrownianMotion(t);
# Specify OU process driven by the Brownian motion
λ = 3.;
σ = 0.1;
ou = OU(λ, σ, bm);
# Generate data from the OU process via exact simulation using Gillespie's algorithm
data = rand(ou, 2.);
In [ ]:
# Plot the simulated OU-BM process
plot(x=t[2:end], y=data, Geom.line)
In [ ]:
# Compute exact MLE estimator of λ using the function based on Brownian motion
# This function does not compute the covariance matrix Σ and its inverse
# BM type is passed to approx_mle_ou_drift
println("Exact MLE estimate of λ via BM: ", exact_mle_ou_drift(bm, data))
# Compute approximate MLE estimator of λ using the function based on FBM
# This function computes the covariance matrix Σ and its inverse
# Since the hurst index h=0.5, FBM coincides with BM and Σ is diagonal
# BM is converted to FBM and resulting FBM is passed to approx_mle_ou_drift
println("Approximate MLE estimate of λ via FBM: ",
approx_mle_ou_drift(convert(FBM, bm), data))
In [ ]:
# Compute approximate MLE estimator of σ using the function based on Brownian motion
# This function does not compute the covariance matrix Σ and its inverse
# BM type is passed to approx_mle_ou_diffusion
println("Approximate MLE estimate of σ via BM: ", approx_mle_ou_diffusion(bm, data))
# Compute approximate MLE estimator of σ using the function based on FBM
# This function computes the covariance matrix Σ and its inverse
# Since the hurst index h=0.5, FBM coincides with BM and Σ is diagonal
# BM is converted to FBM and resulting FBM is passed to approx_mle_ou_diffusion
println("Approximate MLE estimate of σ via FBM: ",
approx_mle_ou_diffusion(convert(FBM, bm), data))
In [ ]:
# Compare exact with approximate MLE estimator of OU-BM
(λ_exact, σ_exact) = exact_mle_ou(bm, data)
println("Exact MLE estimate of (λ, σ) = ", (λ_exact, σ_exact))
(λ_approx, σ_approx) = approx_mle_ou(bm, data)
println("Approximate MLE estimate of (λ, σ) = ", (λ_approx, σ_approx))
In [ ]:
### Simulation 2: comparing exact and approx MLE estimators of OU-BM parameters
### The data are generated by averaging over 10 samples drawn from OU-BM
# Specify Brownian motion
t = 0:1/2^2:5;
bm = BrownianMotion(t);
# Specify OU process driven by the Brownian motion
λ = 3.;
σ = 0.1;
ou = OU(λ, σ, bm);
# Generate data from the OU process via exact simulation using Gillespie's algorithm
# 10 sets of data are generated and their average is stored in data
nsamples = 10;
data = Array(Float64, bm.n-1);
for i = 1:nsamples
data += rand(ou, 2.);
end
data /= nsamples;
In [ ]:
# Plot the average of the 10 simulated OU-BM processes
plot(x=t[2:end], y=data, Geom.line)
In [ ]:
# Compare exact with approximate MLE estimator of OU-BM
(λ_exact, σ_exact) = exact_mle_ou(bm, data)
println("Exact MLE estimate of (λ, σ) = ", (λ_exact, σ_exact))
(λ_approx, σ_approx) = approx_mle_ou(bm, data)
println("Approx MLE estimate of (λ, σ) = ", (λ_approx, σ_approx))
In [ ]:
### Simulation 3: cross-checking the 2 log-likelihood implementations of OU-BM
# Specify Brownian motion
t = 0:1/2^2:5;
bm = BrownianMotion(t);
# Specify OU process driven by the Brownian motion
λ = 3.;
σ = 0.1;
ou = OU(λ, σ, bm);
# Generate data from the OU process via exact simulation using Gillespie's algorithm
data = rand(ou, 2.);
In [ ]:
# Compute approximate MLE estimator of σ of OU-BM given the data
σ_approx = approx_mle_ou_diffusion(bm, data)
# Compute approximate log-likelihood at (λ, σ) = (0.25, σ_approx)
# This function calls uses logpdf of MVNormal
# It is the slower of the 2 ways of computing the log-likelihood
approx_loglik_ou(ou.x, 0.25, σ_approx, data)
In [ ]:
# Compute the components to be passed to alternative log-likelihood implementation
C = autocov(FGN(((bm.t[end]-bm.t[1])/(bm.n-1))^0.5, 0.5), ou.x.n-1);
P = inv(C);
l = [0, data[1:end-1]];
yPy = (data'*P*data)[1];
lPl = (l'*P*l)[1];
yPl = (data'*P*l)[1];
logdetC = logdet(C);
# Compute approximate log-likelihood at (λ, σ) = (0.25, σ_approx)
# This function calls is based on pre-computed data-related functionals
# It is the faster of the 2 ways of computing the log-likelihood
approx_loglik_ou(ou.x, 0.25, approx_mle_ou_diffusion(bm, data), yPy, lPl, yPl, logdetC)
In [ ]:
# Compute the "cached" components to be passed to exact log-likelihood
yy = dot(data, data)
ll = dot(l, l)
yl = dot(data, l)
# Create range of λ values for which the log-likelihoods will be plotted
xrange = 0:1/2^4:6;
xrangelen = length(xrange)-1;
# Compute log-likelihoods for the specified range of λ values
# exact_yvalues holds the values of the exact log-likelihood
# approx_yvalues holds the values of the approximate log-likelihood
exact_yvalues = Array(Float64, xrangelen);
approx_yvalues = Array(Float64, xrangelen);
for i in 1:xrangelen
exact_yvalues[i] = exact_loglik_ou(ou.x, xrange[i+1], σ_approx, yy, ll, yl)
approx_yvalues[i] = approx_loglik_ou(ou.x, xrange[i+1], σ_approx, yPy, lPl, yPl, logdetC)
end
# Compute MLE estimator of λ of OU-BM given data
# This MLE estimate is the maximum of the log-likelihood with respect to λ
exact_mle_ou_drift(bm, data)
println("Exact MLE estimate of λ = ", exact_mle_ou_drift(bm, data))
In [ ]:
plot(x=repeat(collect(xrange[2:end]), outer=[2]),
y=[exact_yvalues, approx_yvalues],
color=repeat(["Exact", "Approximate"], inner=[length(exact_yvalues)]),
Scale.discrete_color_manual("red", "blue"),
Geom.line,
Guide.xlabel("Time points"),
Guide.ylabel("Likelihood"),
Guide.title("Exact vs Approximate Log-likelihood for OU-BM"),
Guide.colorkey("Log-likelihoods"))
In [ ]: