abstract type AbstractGaussianProcess end

Subtypes: GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients. Supertype: Any

Both subtypes represent a Gaussian process as a series expansion in a basis of functions, where the vector of functions is multivariate normally distributed.

See also: GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.

abstract type AbstractModel end

subtype: SDEModel

AbstractModel is a supertype for statistical models.

abstract type AbstactSDE <: Any

subtype: SDE.

A supertype for implementing SDEs.

abstract type AbstractSamplePath <: Any

subtype: SamplePath.

Its subtype implements a continuous samplepath.

struct FaberSchauderExpansionWithGaussianCoefficients{T} <:
        AbstractGaussianProcess where T<:Union{AbstractMvNormal,GaussianVector}
    higestlevel::Int64
    basis::Vector{Function}
    leftboundssupport::Vector{Float64}
    rightboundssupport::Vector{Float64}
    distribution::T
end

Implements a Gaussian Process with Faber-Schauder functions and optimally exploits the sparsity structure of Faber-Schauder functions.

Constructors: FaberSchauderExpansionWithGaussianCoefficients(higestlevel, distribution) FaberSchauderExpansionWithGaussianCoefficients(standarddeviationsperlevel::Vector{Float64})

The user is not allowed to set a basis, which the constructor calculates from the input and is by its very definition a Faber-Schauder basis. The length of distribution should be equal to 2^(higestlevel+1). The second constructor defines a Gaussian process with Faber-Schauder basis with independent coefficients with length(standarddeviationsperlevel)-1 levels (there is also a level zero), so 2^(length(standarddeviationsperlevel)) number of basis functions, where the standard deviation of the coefficients belonging to level k is standarddeviationsperlevel[k+1].

Example

α = 0.5
Π = FaberSchauderExpansionWithGaussianCoefficients([2^(-j*α) for j in 0:5])
f = rand(Π)

struct GaussianProcess{S, T} <: AbstractGaussianProcess where { S<:AbstractVector{U} where U<:Function, T<:Union{AbstractMvNormal,GaussianVector}} basis::S distribution::T end

Implements a Gaussian process defined as an inner product of a Gaussian vector with distribution 'distribution' of type GaussianVector or a subtype of AbstractMvNormal (from the Distributions.jl package) and a vector of functions, the basis of the function space. For the use of a Gaussian process expanded in Faber-Schauder functions up to level j, we recommend using FaberSchauderExpansionWithGaussianCoefficients, as this makes efficiently use of the sparsity structure of Faber-Schauder functions.

See also: FaberSchauderExpansionWithGaussianCoefficients, AbstractGaussianProcess

Example

using Distributions
distribution = MvNormal([k^(-1.0) for k in 1:10])
Π = GaussianProcess([fourier(k) for k in 1:10], distribution)
f = rand(Π)

struct SDE{S, T} <: AbstractSDE where {S<:Function, T<:SDEModel}
  b::S
  model::T
end

Constructors

SDE(b::S, model::T) where {S<:Function, T<:SDEModel}
SDE(model::T, b::S) where {S<:Function, T<:SDEModel}

Implements a stochastic differential equation dXt=b(Xt)dt+σ(Xt)dWt on time interval [0,model.endtime], with Wt a Brownian motion. The drift function b is a function on the real line. model.σ (either a function or a number) gives the diffusion. The begin value X0=model.beginvalue and the sample path is discretised via the Euler-Maruyama scheme with time steps model.Δ.

Example

model = SDEModel(x->2+sin(x), 0.0,10.0,0.001)
sde = SDE(sin, model)

struct SDEModel{T<:Union{Float64, <:Function}} <: AbstractModel
    σ::T
    beginvalue::Float64
    endtime::Float64
    Δ::Float64
end

Constructors

    SDEModel(σ::Q, beginvalue::R, endtime::S, Δ::T) where {Q<:Real, R<:Real,
        S<:Real, T<:Real}
    SDEModel(σ::Q, beginvalue::R, endtime::S, Δ::T) where {Q<:Function, R<:Real,
        S<:Real, T<:Real}

SDEModel implements an SDE model with begin value X_0=beginvalue (at time zero) variance σ, up to end time endtime. The SDE is discretised with precision Δ.

Warning

It is assumed that for every b under consideration, the laws of dXt=b(Xt)dt+σ(Xt)dWt are equivalent.

Examples

SDEModel(1.0, 0.0, 10.0, 0.01)
# Implements the model dX_t=b(X_t)dt+σ(X_t)dW_t, with X_0=0.0, up to time 10,
# with precision 0.01.
#
SDEModel(identity, 0.0, 100.0, 0.1)
# Implements the model dX_t=b(X_t)+X_tdW_t, with X_0=0.0, up to time 100,
# discretised with precision 0.1.

SamplePath(timeinterval::S, samplevalues::T) where
    {S<:AbstractVector{Float64}, T<:AbstractVector{Float64}}

SamplePath implements a continuous samplepath where the timeinterval is not necessarily equally spaced. Sample value samplevalues[k] is the value of the process at time timeinterval[k]. timeinterval is a strictly increasing vector. timeinterval and samplevalues are of equal length.

Performance hint:

Use AbstractRange objects for SamplePath.timeinterval whenever possible.

Constructors:

    SamplePath(timeinterval, samplevalues)
    SamplePath(timeinterval, f::Function) =
        SamplePath(timeinterval, f.(timeinterval))

Examples

t = 0.0:0.01:2.0
X = SamplePath(t, sinpi)

t = 0.0:0.01:1.0
v = map(x->floor(10*x), t)
X = SamplePath(t, v)

Base.eltype(::Type{SamplePath})

Outputs element type of SamplePath, which is Float64.

Base.firstindex(X::SamplePath)

The first index of the SamplePath is 1.

Base.getindex(X::SamplePath, i::Int)

Returns the i-th observation of the samplepath, which is X.samplevalues[i].

Base.iterate(X::SamplePath, state=1)

Iterate over the sample values.

Example

X = SamplePath([1.0, 2.0, 3.0], x->x^2)
for value in X
    println(value)
end
# Or in the reverse
for value in Iterators.Reverse(X)
    println(value)
end 

Base.lastindex(X::SamplePath)

The last index of the SamplePath is length(X).

length(Π::AbstractGaussianProcess)
length(Π::FaberSchauderExpansionWithGaussianCoefficients)

length returns the number of basis functions, which is equal to the length of the distribution of the coefficients, which is also equal to the number of coefficients. For a FaberSchauderExpansionWithGaussianCoefficients object, the length is equal to 2^(higestlevel+1).

Example

Π = GaussianProcess([sin, cos], MvNormal([1.,1.]))
length(Π)
#
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^(-0.5*j) for j in 0:3])
length(Π) # == 2^4 == 16.

Base.length(X::SamplePath)

Returns the length of the vector timeinterval == length vector samplevalues.

Examples

X = SamplePath([0.,1.,2.], [3.,5., -1.])
length(X) # == 3

X = SamplePath(0.0:0.1:2π, sin)
length(X) == length(0.0:0.1:2π)

rand(Π::AbstractGaussianProcess)

rand returns a random function, where the coefficients have distribution Π.distribution and the basis functions are defined in Π.basis.

Example

using Distributions
distribution = MvNormal([k^(-1.0) for k in 1:100])
Π = GaussianProcess([fourier(k) for  k in 1:100], distribution)
f = rand(Π)

rand(sde::SDE)

Returns a sample path from the SDE sde, from time 0.0 to time sde.model.endtime, discretised with precision sde.model.Δ with the help of the Euler-Maruyama scheme.

Examples

model = SDEModel(1.0,0.0,10.0,0.01)
sde = SDE(sin, model)
X = rand(sde)

Base.step(X::SamplePath{<:AbstractRange})

Returns the step size of SamplePath.timeinterval.

Example

t = 0.0:0.01:2.0
X = SamplePath(t, sinpi)
step(X) == step(t)

calculateboundssupport(higestlevel::Int)

Internal function, not exported.

Returns a tuple of two Float64 vectors of length 2^(higestlevel+1) where the i-th element of the first vector is the left bound of the support of the i-th Faber-Schauder basis function and the i-th element of the second vector is the right bound of the support of i-th Faber-Schauder basis function.

calculatedependentfaberschauderfunctions(higestlevel::Int64)

Internal function, not exported.

Calculates which Faber-Schauder functions (ψi,ψj) have essential overlapping support, i ≤ j. All other combinations (ψi, ψj), i ≤ j, have essential nonoverlapping support and hence int0^T ψi(Xt)ψj(X_t)dt=0, which we use in calculategirsanovmatrix.

Returns a triple (lengthvectors, rowindices, columnindices) where

lengthvectors == length(rowindices) == length(columnindices),

and for i in 1:lengthvectors (ψrowindices[i],ψcolumnindices[i]) have essentially overlapping support. For all i, rowindices[i] ≤ columnindices[i].

calculategirsanovmatrix(sizesquarematrix, timeinterval, ψXt, σXt)
calculategirsanovmatrix(Π, samplevalueindices, timeinterval, ψXt, σXt)

Internal function, not exported!

Calculates the Girsanov matrix. In the case of a FaberSchauderExpansionWithGaussianCoefficients-object it will make use of the sparsity structure of the Faber-Schauder basis.

calculategirsanovmatrixelement(ψ1Xt, ψ2Xt, σ, Δt)
calculategirsanovmatrixelement(samplevalueindices1, samplevalueindices2, ψ1Xt, ψ2Xt, σ, Δt)

Internal function, not exported!

Where ψ1Xt, ψ2Xt are arrays, and σ and Δt are either arrays or numbers.

Let ψa and ψb denote two basis elements. calculategirsanovmatrixelement calculates the (a,b) element of the Girsanov matrix defined by int0^T ψa(Xt)ψb(Xt)/(σ^2(Xt)) dt, where ψa is given by ψ1Xt (already evaluated in ψa) and ψb by ψ2Xt (already evaluated in ψb).

samplevalueindices1[i] == false should correspond to Ψ1Xt[i] == 0.0. Similar for samplevalueindices2 and Ψ2Xt.

calculategirsanovvector(lengthvector, samplevalues, ψXt, σXt)
calculategirsanovvector(lengthvector, samplevalueindices, samplevalues, ψXt, σXt)

Internal function, not exported!

Calculates the Girsanov vector.

calculategirsanovvectorelement(ΔXt, ψXt, σXt)
calculategirsanovvectorelement(ΔXt, ψXt, σXt::Number)

calculategirsanovvectorelement(samplevalueindices, ΔXt, ψXt, σ::Number)

Internal function, used by calculategirsanovvector, not exported!

Calculates the kth Girsanov vector element int0^T ψk(Xt)/(σ^2(Xt)) dXt with the Euler-Maruyama method where ψk is the kth basis function, and T is the end time.

calculatenextsamplevalue(prevXval, b, model::SDEModel{<:Number}, BMincrement)
calculatenextsamplevalue(prevXval, b, model, BMincrement)

Internal function, not exported.

Given the previous samplepath value prevXval, it calculates the next sample value with the Euler-Maruyama scheme.

calculateposterior(Π, X, model::SDEModel)
calculateposterior(Π, X, σ)
calculateposterior(Π::FaberSchauderExpansionWithGaussianCoefficients, X, σ)

Calculates the posterior distribution Π(⋅∣X) and returns a FaberSchauderExpansionWithGaussianCoefficients object when Π is a FaberSchauderExpansionWithGaussianCoefficients-object. Otherwise, it returns a GaussianProcess-object. Uses model to determine σ. When Π is a FaberSchauderExpansionWithGaussianCoefficients-object, then calculateposterior makes use of the sparsity structure of the Faber-Schauder basis.

Examples

##Example with Faber-Schauder expansion.

model = SDEModel(1.0, 0.0, 10000.0, 0.1)
sde = SDE(x->sinpi(2*x), model)
X = rand(sde)
Π = FaberSchauderExpansionWithGaussianCoefficients([2.0^j for j in 0:5])
postΠ = calculateposterior(Π, X, model)

##Example with Fourier expansion.

using Distributions
α = 0.5
model = SDEModel(1.0, 0.0, 10000.0, 0.1)
sde = SDE(x->sinpi(2*x), model)
X = rand(sde)
distribution = MvNormalCanon([k^(α+0.5) for k in 1:50])
Π = GaussianProcess([fourier(k) for k in 1:50], distribution)
postΠ = calculateposterior(Π, X, model)

calculateΔt(timeinterval::AbstractRange) = step(timeinterval)
calculateΔt(timeinterval) = timeinterval[2:end] - timeinterval[1:end-1]

Internal function, not exported!

Calculates Δt. It returns a number when timeinterval is a range object and an array of increments of the timeinterval vector otherwise.

calculateσXt(σ::Number, v) = σ
calculateσXt(σ, v) = σ.(v)

Internal function, not exported!.

Returns σ when it is a number, and if σ is a function, it evaluates v in σ.

createFaberSchauderBasisUpToLevelHigestLevel(higestlevel::Int)

Internal function, not exported.

Creates Faber-Schauder basis up to level higestlevel. Returns a Vector{Function} type of length 2^(higestlevel+1). The first element of the vector is faberschauderone, and the 2^j+k element is faberschauder(j,k), with 0≤j≤higestlevel and 1≤k≤2^j.

createvectorofstandarddeviationsfromstandarddeviationsperlevel(standarddeviationsperlevel::AbstractArray{Float64})

Internal function, not exported!

Calculated a vector of standard deviations for all basis functions, where basis function i gets standard deviation standarddeviationsperlevel[ℓ+1] if ℓ is the level of i (we start at level 0).

faberschauder(j::Int, k::Int)

faberschauder implements the k-th Faber-Schauder function of level j. Here, j≥0 and 1≤k≤2^j. It is a one-periodic function and defined on [0,1] by 2^(j+1)(x-(k-1)2^(-j)) on (k-1)2^(-j)≤x≤(k-1/2)2^(-j) and 1 - 2^(j+1)(x-(k-1/2)2^(-j)) on [(k-1/2)2^(-j), k2^(-j)] and zero outside these intervals.

See also: faberschauderone.

#Warning

Note the difference between faberschauderone and faberschauder. The first is a function that takes a Float64 and returns a Float64, the second takes (j,k) and returns an anonymous function that takes a Float64 and returns a Float64.

#Example with Plots

using Plots
J=2
x=0.0:0.001:1.0
p=plot()
for j in 0:J
  for k in 1:2^j
      y = faberschauder(j,k).(x)
      plot!(p,x,y)
  end
end
display(p)

faberschauderone(x)

Implements the first Faber-Schauder function defined by 1-2x for 0≤x≤1/2 and -1+2x for 0.5≤x≤1, and is 1-periodically extended to all x∈R.

#See also: faberschauder

#Warning

Note the difference between faberschauderone and faberschauder. The first is a function that takes a Float64 and returns a Float64, the second takes (j,k) and returns an anonymous function that takes a Float64 and returns a Float64.

#Examples

x=-2.0:0.001:2.0
y=faberschauderone.(x)

fourier(k::Int)

fourier implements the Fourier basis of functions ϕk, defined by ϕ0≡1, and if k≥1 and k is odd, by ϕk(x)=sqrt(2)sin((k+1)πx) and if k≥1 and even, by ϕk(x)=sqrt(2)cos(πkx). This function is defined for k≥0.

Examples

x = -1.0:0.01:2.0
y = fourier(3).(x)

isstrictlyincreasing(x::AbstractVector{T})::Bool where T <: Number
isstrictlyincreasing(x::AbstractRange{T})::Bool where T <: Number

Internal function, not exported.

Tests whether a vector of numbers is strictly increasing.

sumoffunctions(vectoroffunctions::AbstractArray{<:Function},
    vectorofscalars::AbstractArray{<:Number})

Internal function, not exported.

sumoffunctions calculates the 'inner product' of the function vector and the scalar vector. In other words, the sum of the functions weighted by vectorofscalars. It returns a function.

mean(Π::AbstractGaussianProcess)

mean calculates the mean of the Gaussian process. It returns a function.

#Examples

using Distributions
distribution = MvNormal([k^(-1.0) for k in 1:100])
Π = GaussianProcess([fourier(k) for  k in 1:100], distribution)
mean(Π)