Getting Started
Installation
WhittleLikelihoodInference.jl can be installed by running the following command:
] add https://github.com/JakeGrainger/WhittleLikelihoodInference.jlThe package can then be used by running
using WhittleLikelihoodInferenceA basic example (univariate)
For a simple univariate example, consider a stationary OU process. The stationary OU process has parameters $\theta$ and $\sigma$ with autocovariance
\[c(\tau) = \frac{\sigma^2}{\theta}\exp\{-\theta|\tau|\}\]
and spectral density function
\[f(\omega) = \frac{\sigma^2}{\pi(\theta^2+\omega^2)}.\]
The OU process is a continuous time Gaussian process, which we are interested in fitting to a sampled process which has finite length. We can simulate such a process with length $n$ and sampling interval $\Delta$ with the following code:
σ = 2
θ = 0.5
n = 1000
Δ = 0.5
nreps = 1
ts = simulate_gp(OU(σ,θ),n,Δ,nreps)[1]WhittleLikelihoodInference.TimeSeries{Float64, 1}([2.4532283525716148, 0.310710140980738, -0.6353185948906841, -2.0967292921396816, -0.09932113629970596, 3.8469424085764734, 3.941314442967764, 2.587388202543295, 2.9062974448223597, 1.346198890593091 … -0.19740090802332655, 0.611584915089268, -1.5652354777080935, -1.7119898043444108, -2.5803699473209, -1.2231151156598878, -3.3581289784651616, -1.3988162008015068, 0.18276780465571485, -0.0831926047571335], 0.5)In this case, we only wanted one series, but simulate_gp will simulate a vector of TimeSeries of length nreps. Therefore, we pull out the first index of this array as we only want one series.
Alternatively, simulation can be performed by running the following:
import WhittleLikelihoodInference: GaussianProcess
X = GaussianProcess(OU(σ,θ),n,Δ)
ts = rand(X)WhittleLikelihoodInference.TimeSeries{Float64, 1}([2.4532283525716148, 0.310710140980738, -0.6353185948906841, -2.0967292921396816, -0.09932113629970596, 3.8469424085764734, 3.941314442967764, 2.587388202543295, 2.9062974448223597, 1.346198890593091 … -0.19740090802332655, 0.611584915089268, -1.5652354777080935, -1.7119898043444108, -2.5803699473209, -1.2231151156598878, -3.3581289784651616, -1.3988162008015068, 0.18276780465571485, -0.0831926047571335], 0.5)To perform Whittle likelihood inference we need to first create the objective function. This is done with the following:
whittle_objective = WhittleLikelihood(OU,ts)Whittle likelihood for the OU model.The resulting objective function is set up to work with Optim so we can use
julia> whittle_objective([2, 0.5])-372.04805671473946julia> F,G,H = 0,zeros(2),zeros(2,2)(0, [0.0, 0.0], [0.0 0.0; 0.0 0.0])julia> whittle_objective(F,G,H,[2, 0.5])-372.04805671473946julia> G2-element Vector{Float64}: -31.387533331731234 16.992769821055777julia> H2×2 Matrix{Float64}: 1031.39 -477.041 -477.041 993.791
We can optimise this function with Optim by using
using Optim
constraints = Optim.TwiceDifferentiableConstraints(zeros(2),fill(Inf,2))
obj = TwiceDifferentiable(Optim.only_fgh!(whittle_objective),ones(2))
res = optimize(obj, constraints, [1.5, 0.2], IPNewton()) * Status: success
* Candidate solution
Final objective value: -3.725345e+02
* Found with
Algorithm: Interior Point Newton
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 1.06e-06 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 31
f(x) calls: 190
∇f(x) calls: 190
The full example is:
using WhittleLikelihoodInference, Optim
import WhittleLikelihoodInference: GaussianProcess
σ = 2
θ = 0.5
n = 1000
Δ = 0.5
ts = rand(GaussianProcess(OU(σ,θ),n,Δ))
whittle_objective = WhittleLikelihood(OU,ts)
constraints = Optim.TwiceDifferentiableConstraints(zeros(2),fill(Inf,2))
obj = TwiceDifferentiable(Optim.only_fgh!(whittle_objective),ones(2))
res = optimize(obj, constraints, [1.5, 0.2], IPNewton()) * Status: success
* Candidate solution
Final objective value: -3.725345e+02
* Found with
Algorithm: Interior Point Newton
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 1.06e-06 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 31
f(x) calls: 190
∇f(x) calls: 190
A basic example (multivariate)
An analagous example for a multivariate process is a bivariate correlated OU process. Consider two independent OU processes $X_1$ and $X_2$ with the same $\sigma$ and $\theta$ parameters. Then define $Y_1 = X_1$ and $Y_2 = \rho X_1 + \sqrt{1-\rho^2}X_2$. We have that
\[\mathbb{E}[Y_1(\tau)Y_1(0)] = \mathbb{E}[X_1(\tau)X_1(0)]\]
and
\[\begin{align*} \mathbb{E}[Y_2(\tau)Y_2(0)] &= \rho^2\mathbb{E}[X_1(\tau)X_1(0)] + (1-\rho^2)\mathbb{E}[X_2(\tau)X_2(0)]\\ &= \mathbb{E}[X_1(\tau)X_1(0)] \end{align*}\]
and
\[\begin{align*} \mathbb{E}[Y_2(\tau)Y_1(0)] &= \rho\mathbb{E}[X_1(\tau)X_1(0)] \end{align*}.\]
So we have taken two OU processes and correlated them with correlation parameter $\rho$. We can fit the Whittle likelihood to a simulated example with the below code.
using WhittleLikelihoodInference, Optim
import WhittleLikelihoodInference: GaussianProcess
σ = 2
θ = 0.5
ρ = 0.7
n = 1000
Δ = 0.5
ts = rand(GaussianProcess(CorrelatedOU(σ,θ,ρ),n,Δ))
whittle_objective = WhittleLikelihood(CorrelatedOU,ts)
constraints = Optim.TwiceDifferentiableConstraints(zeros(3),[Inf,Inf,1])
obj = TwiceDifferentiable(Optim.only_fgh!(whittle_objective),ones(3))
res = optimize(obj, constraints, [1.5, 0.2, 0.5], IPNewton()) * Status: success
* Candidate solution
Final objective value: -1.420219e+03
* Found with
Algorithm: Interior Point Newton
* Convergence measures
|x - x'| = 0.00e+00 ≤ 0.0e+00
|x - x'|/|x'| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00
|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
|g(x)| = 5.18e-05 ≰ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 34
f(x) calls: 213
∇f(x) calls: 213