The YUIMA Project is an open source academic project aimed at developing a complete environment for estimation and simulation of Stochastic Differential Equations and other Stochastic Processes via the R package called yuima
and its Graphical User Interface yuimaGUI
.
The YUIMA Object
The main object is the yuima
object which allows to describe the model in a mathematically sound way. Then the data and the sampling structure can be included as well for estimation and simulation purposes.
Model
The ‘setModel’ function
The setModel()
function defines a stochastic differential equation with or without jumps of the following form:
\[ dX_t = a(t,X_t, \alpha)dt + b(t,X_t,\beta)dW_t^H + c(t,X_t,\gamma)dZ_t \] where
- \(a(t,X_t,\alpha)\) is the drift term. Described by the
drift
argument - \(b(t,X_t,\beta)\) is the diffusion term. Described by the
diffusion
argument - \(c(t,X_t,\gamma)\) is the jump term. Described by the
jump.coeff
argument - \(H\) is the Hurst coefficient. Described by the
hurst
argument - \(Z_t\) is the Levy noise. Described by the
measure.type
andmeasure
arguments
Deterministic Model
\[ dU_t = \sin(\alpha t) dt \]
setModel(drift = "sin(alpha*t)", # the drift term
solve.variable = "u", # the solve variable
time.variable = "t") # the time variable
Geometric Brownian Motion
\[ dX_t = \mu X_t \; dt + \sigma X_t \; dW_t \]
setModel(drift = "mu*x", # the drift term
diffusion = "sigma*x", # the diffusion term
solve.variable = "x") # the solve variable
CKLS Model
\[ dX_t = (\theta_1+\theta_2 X_t) \; dt + \theta_3 X_t^{\theta_4} \; dW_t \]
setModel(drift = "theta1+theta2*x", # the drift term
diffusion = "theta3*x^theta4", # the diffusion term
solve.variable = "x") # the solve variable
2-Dimensional Diffusion with 3 Noises
\[ \begin{cases} dX_t^1 = -3X_t^1 \; dt + dW_t^1 + X_t^2 dW_t^3 \\ dX_t^2 = -(X_t^1+2X_t^2) \; dt + X_t^1 dW_t^1 + 3 dW_t^2 \end{cases} \]
setModel(drift = c("-3*x1","-x1-2*x2"), # the drift vector
diffusion = matrix(c("1","x1","0","3","x2","0"), 2, 3), # the diffusion matrix
solve.variable = c("x1","x2")) # the solve variables
Fractional Ornstein-Uhlenbeck
\[ dX_t = -\theta X_t \; dt + \sigma \; dW_t^H \]
setModel(drift = "-theta*x", # the drift term
diffusion="sigma", # the diffusion term
hurst = NA, # the hurst coefficient
solve.variable = "x") # the solve variable
Jump Process with Compound Poisson Measure
\[ dX_t = -\theta X_t dt + \sigma dW_t + dZ_t \]
setModel(drift = "-theta*x", # the drift term
diffusion="sigma", # the diffusion term
jump.coeff = "1", # the jump term
measure.type = "CP", # the measure type
measure = list( # the measure
intensity = "lambda", # constant intensity
df = "dnorm(z, mu_jump, sigma_jump)" # jump density function
),
solve.variable = "x") # the solve variable
The ‘setPoisson’ Function
Defines a generic Compound Poisson model.
Compound Poisson with constant intensity and Gaussian jumps
\[ X_t = X_0+\sum_{i=0}^{N_t} Y_i \; : \;\;\; N_t \sim Poi\Bigl(\int_0^t \lambda(t)dt\Bigl) , \;\;\;\; Y_i \sim N(\mu_{jump}, \; \sigma_{jump}) \\ \lambda(t)=\lambda \]
setPoisson(intensity = "lambda", # the intensity function
df = "dnorm(z, mean = mu_jump, sd = sigma_jump)", # the density function
solve.variable = "x") # the solve variable
Compound Poisson with exponentially decaying intensity and Student-t jumps
\[ X_t = X_0+\sum_{i=0}^{N_t} Y_i \; : \;\;\; N_t \sim Poi\Bigl(\int_0^t \lambda(t)dt\Bigl) , \;\;\;\; Y_i \sim t( \nu_{jump}, \; \mu_{jump} ) \\ \lambda(t)=\alpha \; e^{-\beta t} \]
The ‘setCarma’ Function
Defines a generic Continuous ARMA model.
Continuous ARMA(3,1) process driven by a Brownian Motion \[ CARMA(3,1) \]
Continuous ARMA(3,1) process driven by a Compound Poisson with Gaussian jumps \[ CARMA(3,1) \]
The ‘setCogarch’ Function
Defines a generic Continuous GARCH model.
Continuous COGARCH(1,1) process driven by a Compound Poisson with Gaussian jumps \[ COGARCH(1,1) \]
Data
The setData()
function prepares the data for model estimation. The delta
argument describes the time increment between observations. If we have monthly data and want to measure time in years, then delta
should be \(1/12\). If we have daily data and want to measure time in months, then delta
should be \(1/30\). If we have financial daily data and want to measure time in years, then delta
should be \(1/252\), since 252 is the average number of trading days in one year. In general, if we want to measure time in unit \(T\), delta
should be 1 over the average number of observations in a period \(T\). The unit of measure of time affects the estimated value of the model parameters.
The following example downloads and sets some financial data (see tutorial on Data Acquisition in R).
# Install the quantmod package if needed:
# install.packages('quantmod')
# load quantmod
require(quantmod)
# download Facebook quotes
fb <- getSymbols(Symbols = 'FB', src = 'yahoo', auto.assign = FALSE)
# setData with time in years -> delta = 1/252
# (there are 252 observations in 1 year)
setData(fb$FB.Close, delta = 1/252, t0 = 0)
##
##
## Number of original time series: 1
## length = 2016, time range [2012-05-18 ; 2020-05-22]
##
## Number of zoo time series: 1
## length time.min time.max delta
## FB.Close 2016 0 7.996 0.003968254
Sampling
The setSampling()
function describes the simulation grid. If delta
is not specified, it is calculated as (Terminal-Initial)/n
. If delta
is specified, the Terminal
is adjusted to be equal to Initial+n*delta
.
Simulation
Simulation of a generic model is perfomed with the simulate()
function.
Example Solve an Ordinary Differential Equation
# model: ordinary differential equation
model <- setModel(drift = 'sin(t)*t', solve.variable = 'x', time.variable = 't')
# simulation scheme
sampling <- setSampling(Initial = 0, Terminal = 10, n = 1000)
# yuima object
yuima <- setYuima(model = model, sampling = sampling)
# simulation
sim <- simulate(yuima)
# plot
plot(sim)
Example Simulate one trajectory of a jump diffusion model
# model: jump diffusion
model <- setModel(drift = "-theta*x",
diffusion="sigma",
jump.coeff = "1",
measure.type = "CP",
measure = list(
intensity = "lambda",
df = "dnorm(z, mu_jump, sigma_jump)"
),
solve.variable = "x")
# simulation scheme
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)
# yuima object
yuima <- setYuima(model = model, sampling = sampling)
# simulation
sim <- simulate(yuima, # the yuima object
xinit = 1, # the initial value
true.parameter = list( # specify the parameters:
theta = 1, # value for the 'theta' parameter
sigma = 1, # value for the 'sigma' parameter
lambda = 10, # value for the 'lambda' parameter
mu_jump = 0, # value for the 'mu_jump' parameter
sigma_jump = 2 # value for the 'sigma_jump' parameter
))
# plot
plot(sim)
Estimation
The qmle()
function calculates the quasi-likelihood and estimate of the parameters of the stochastic differential equation by the maximum likelihood method or least squares estimator of the drift parameter.
Example Simulate a Geometric Brownian Motion and estimate its parameters
# model: geometric brownian motion
model <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')
# simulation scheme
sampling <- setSampling(Initial = 0, Terminal = 1, n = 1000)
# yuima object
yuima <- setYuima(model = model, sampling = sampling)
# simulation
sim <- simulate(yuima, true.parameter = list(mu = 1.3, sigma = 0.25), xinit = 100)
# estimation
estimation <- qmle(sim, # the yuima object
start = list(mu = 0, sigma = 1), # starting values for optimization
lower = list(sigma = 0)) # lower bounds
# estimates and standard errors
summary(estimation)
## Quasi-Maximum likelihood estimation
##
## Call:
## qmle(yuima = sim, start = list(mu = 0, sigma = 1), lower = list(sigma = 0))
##
## Coefficients:
## Estimate Std. Error
## sigma 0.2482216 0.005630911
## mu 1.1125000 0.248221634
##
## -2 log L: 3150.704
Example Estimate the yearly volatility (\(\sigma\) in the Geometric Brownian Motion) of Google stock quotes
# Install the quantmod package if needed:
# install.packages('quantmod')
# load quantmod
require(quantmod)
# download Google quotes
goog <- getSymbols(Symbols = 'GOOG', src = 'yahoo', auto.assign = FALSE)
# setData with time in years -> delta = 1/252
# (there are 252 observations in 1 year)
data <- setData(goog$GOOG.Close, delta = 1/252, t0 = 0)
# model: geometric brownian motion
model <- setModel(drift = 'mu*x', diffusion = 'sigma*x', solve.variable = 'x')
# yuima object
yuima <- setYuima(model = model, data = data)
# estimation
estimation <- qmle(yuima, # the yuima object
start = list(mu = 0, sigma = 0.5), # starting values for optimization
lower = list(sigma = 0)) # lower bounds
# estimates and standard errors
summary(estimation)
## Quasi-Maximum likelihood estimation
##
## Call:
## qmle(yuima = yuima, start = list(mu = 0, sigma = 0.5), lower = list(sigma = 0))
##
## Coefficients:
## Estimate Std. Error
## sigma 0.2939155 0.003581592
## mu 0.1781250 0.080372571
##
## -2 log L: 24191.29
yuimaGUI
The yuimaGUI package provides a user-friendly interface for yuima. It simplifies tasks such as estimation and simulation of stochastic processes, including additional tools related to quantitative finance such as data retrieval of stock prices and economic indicators, time series clustering, change point analysis, lead-lag estimation.
The yuimaGUI is available online for free, but it is strongly recommended to install the application via the R package on your local machine for better performance and less downtime.