We use the historic data to estimate the parameters for Ornstein Uhlenbeck process of crop prices in Alberta. And then we use the estimated parameters for Monte Carlo simulation.
The Ornstein-Uhlenbeck or Vasicek process is the unique solution to the following stochastic differential equation:(Stochastic Differential Equation, 2008, p44.)
Ornstein-Uhlenbeck model \[ dX_t = - \theta_2 X_t dt +\theta_3dW_t, \]
Vasicek modified it to
\[ dX_t = (\theta_1 - \theta_2X_t)dt +\theta_3dW_t, X_0 = x_0 \]
\( (\theta_1 - \theta_2X_t)dt \) is deterministic part; \( \theta_3dW_t \) is stochastic part.
\( dW_t \) is the Brownian motion, which follows random normal distributioin \( N(0,t) \).
\( (\theta_1 - \theta_2X_t) \) is drift; \( \theta_3 \) is diffusion.
In finance, the model more often is written as:
\[ dS_t = \theta(\mu - S_t)dt + \sigma dW_t \]
where \( \sigma \) is interpreted as the instantaneous volatility , \( \sigma^2/(2\theta) \) is the long term variance; \( \mu \) is the long-run equilibrium value of the process, and \( \theta \) is the speed of reversion.
\( \theta \) increases the speed at which the system stabilizes around the long term mean \( \mu \).
\( \sigma \) increases the amount of randomness entering the system.
For \( \theta_2>0 \), \( X_t~ \) are assumed iid \[ N(\frac{\theta_1}{\theta_2},\frac{\theta_3^{2}}{2\theta_2}) \].
\[ \frac{\theta_1}{\theta_2} = \mu ~~\mbox{and}~~ \theta_2 = \theta ~~\mbox{and}~~ \theta_3 = \sigma \].
For any \( t>=0 \), the density of the distribution of \( X_t \) given \( X_0 = x_0 \), with mean and variance respectively:
\[ m(t,x) = E_\theta(X_t ~|~ X_0 = x_0) = \frac{\theta_1}{\theta_2} + (x_0-\frac{\theta_1}{\theta_2})e^{-\theta_2t} \]
and
\[ v(t,x) = Var_\theta(X_t ~|~ X_0 = x_0) =\frac{\theta_3^2 (1- e^{-2 \theta_2 t} )}{2 \theta_2} \]
We got monthly price data of crops in Alberta from Statistics Canadahttp://www5.statcan.gc.ca/subject-sujet/result-resultat?pid=920&id=2024&lang=eng&type=ARRAY&sortType=3&pageNum=0.
setwd("E:/Dropbox/book/economics/485/projects/pricemodel/oumodel")
crop.price <- read.csv("abcropprice.csv", header = T, sep = ",")
# set the date format
crop.price[, 1] <- as.Date(crop.price[, 1], format = "%d/%m/%Y")
head(crop.price)
## Date Wheat.excluding.durum Oats Barley Canola Dry.peas
## 1 1985-01-01 166.0 114.6 127.4 342.0 NA
## 2 1985-02-01 167.2 116.2 127.4 347.3 NA
## 3 1985-03-01 166.9 118.4 130.3 350.0 NA
## 4 1985-04-01 164.8 115.7 127.4 363.8 NA
## 5 1985-05-01 166.8 105.7 125.6 354.6 NA
## 6 1985-06-01 166.6 102.2 121.4 349.6 NA
For the Ornstein Uhlenbeck Model, long term data has more breaks or jumps, so we only use the short term data which is from stable period and capture the current characteristics of the stochastic process. We decide to use the data after 2007.
price <- crop.price[crop.price$Date > "2006-12-01", ]
head(price)
## Date Wheat.excluding.durum Oats Barley Canola Dry.peas
## 265 2007-01-01 169.9 149.0 134.5 339.1 181
## 266 2007-02-01 170.7 152.6 142.7 350.1 189
## 267 2007-03-01 174.5 163.4 144.8 353.0 205
## 268 2007-04-01 167.2 165.3 147.2 354.0 230
## 269 2007-05-01 165.2 165.3 154.4 362.6 239
## 270 2007-06-01 167.4 167.4 160.8 370.8 249
There are many ways to calculate the parameters. We tried two ways. First, we use package “sde.sim” (Stochastic Differential Equation, 2008);
Recap \[ \frac{\theta_1}{\theta_2} = \mu ~~\mbox{and}~~ \theta_2 = \theta ~~\mbox{and}~~ \theta_3 = \sigma \].
Get the necessary packeages and functions
# install.packages("stats4")
# install.packages("sde")
require ( stats4 )
require ( sde )
# Two functions from page 114
# Functions dcOU evaluates the conditional density of the process. It generates the density of the x which from cerntain distribution
dcOU<- function (x, t, x0 , theta , log = FALSE ){
Ex <- theta [1] / theta [2]+( x0 - theta [1] / theta [2]) * exp (- theta [2] *t)
Vx <- theta [3]^2 *(1- exp (-2* theta [2] *t))/(2* theta [2])
dnorm (x, mean =Ex , sd = sqrt (Vx), log = log )
}
# Function OU.lik generates the log likelihood function of X for the MLE process.
OU.lik <- function ( theta1 , theta2 , theta3 ){
n <- length (X)
# deltat is the time interval between observations
dt <- deltat (X)-sum (dcOU (X [2: n], dt , X [1:(n -1)] , c( theta1 , theta2 , theta3 ), log = TRUE ))
}
# The function OU.lik needs as input the three parameters and assumes that sample observations of the process are available in the current R workspace in the ts object X.
Mle is very sensitive for the start value. We choose theta1 =20, theta2 =0.1 , theta3 =1. We cannot get it converg. It seems not working. (The result value which the process converges to might not be the unique solution. )
wheat <- as.vector(price[,2]) # convert to vector
head(wheat)
X <- ts(data=wheat)
mle(OU.lik , start = list ( theta1 =20, theta2= 0.1 , theta3 =10) , method ="L-BFGS-B", lower =c(-Inf ,0 ,0)) -> fitwheat
summary ( fitwheat )
coef(fitwheat)[[1]]/coef(fitwheat)[[2]]
Second, we manually use the function from “Calibrating the Ornstein-Uhlenbeck (Vasicek) model”http://www.sitmo.com/article/calibrating-the-ornstein-uhlenbeck-model/
Function “ouFit.ML” is saved in file “oufit.R”.
# function for Calibration using Maximum Likelihood estimates
ouFit.ML = function(spread) {
n = length(spread)
delta = 1 # delta
Sx = sum(spread[1:n - 1])
Sy = sum(spread[2:n])
Sxx = sum((spread[1:n - 1])^2)
Syy = sum((spread[2:n])^2)
Sxy = sum((spread[1:n - 1]) * (spread[2:n]))
mu = (Sy * Sxx - Sx * Sxy)/((n - 1) * (Sxx - Sxy) - (Sx^2 - Sx * Sy))
theta = -log((Sxy - mu * Sx - mu * Sy + (n - 1) * mu^2)/(Sxx - 2 * mu *
Sx + (n - 1) * mu^2))/delta
a = exp(-theta * delta)
sigmah2 = (Syy - 2 * a * Sxy + a^2 * Sxx - 2 * mu * (1 - a) * (Sy - a *
Sx) + (n - 1) * mu^2 * (1 - a)^2)/(n - 1)
sigma = sqrt((sigmah2) * 2 * theta/(1 - a^2))
theta = list(theta = theta, mu = mu, sigma = sigma, sigmah2 = sigmah2)
return(theta)
}
source("oufit.R")
wheat <- as.vector(price[, 2]) # convert to vector
head(wheat)
## [1] 169.9 170.7 174.5 167.2 165.2 167.4
plot(wheat, type = "l")
ouFit.ML(wheat)
## $theta
## [1] 0.1198
##
## $mu
## [1] 249.9
##
## $sigma
## [1] 21.99
##
## $sigmah2
## [1] 430.1
oats <- as.vector(na.omit(price[, 3])) # get rid of n/a and convert to vector
head(oats)
## [1] 149.0 152.6 163.4 165.3 165.3 167.4
plot(oats, type = "l")
ouFit.ML(oats)
## $theta
## [1] 0.08256
##
## $mu
## [1] 186.2
##
## $sigma
## [1] 9.471
##
## $sigmah2
## [1] 82.68
barley <- as.vector(na.omit(price[, 4])) # get rid of n/a and convert to vector
head(barley)
## [1] 134.5 142.7 144.8 147.2 154.4 160.8
plot(barley, type = "l")
ouFit.ML(barley)
## $theta
## [1] 0.04971
##
## $mu
## [1] 188.8
##
## $sigma
## [1] 11.29
##
## $sigmah2
## [1] 121.3
canola <- as.vector(na.omit(price[, 5])) # get rid of n/a and convert to vector
head(canola)
## [1] 339.1 350.1 353.0 354.0 362.6 370.8
plot(canola, type = "l")
ouFit.ML(canola)
## $theta
## [1] 0.04689
##
## $mu
## [1] 507.7
##
## $sigma
## [1] 18.46
##
## $sigmah2
## [1] 325.5
dry.peas <- as.vector(na.omit(price[, 6])) # get rid of n/a and convert to vector
head(dry.peas)
## [1] 181 189 205 230 239 249
plot(dry.peas, type = "l")
ouFit.ML(dry.peas)
## $theta
## [1] 0.05435
##
## $mu
## [1] 279.5
##
## $sigma
## [1] 14.08
##
## $sigmah2
## [1] 187.7
There are two packages we can use to generate Vasicek process. We try both and get different results. It may be because both use the random method to simulation.
We find the “yuima” is better than “sde.sim”
The simulation follows the method mentioned on “1 4 Vasicek Model ”https://www.youtube.com/watch?v=5BpOYPNxWsA&index=4&list=PLRj8HuwfWZEsXW2pzAwAWYC8EZbD2ieOq
Recall “yuima” package uses this notation:
\[ dS_t = \theta(\mu - S_t)dt + \sigma dW_t \]
The deltat \( \dt \) is set \( =1/n=1/100 \) by default. We can change it by using \( grid=setSampling(Terminal=1,n=1000) \).
One simulation, different time different result. And time period still is 1000 (“yuimu” default: delta = 1/100, we change it to 1/1000).
require(yuima)
## Loading required package: yuima
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: stats4
## Loading required package: expm
## Loading required package: Matrix
##
## Attaching package: 'expm'
##
## The following object is masked from 'package:Matrix':
##
## expm
##
## ############################################
## This is YUIMA Project package.
## Check for the latest development version at:
## http://R-Forge.R-Project.org/projects/yuima
## ############################################
##
## Attaching package: 'yuima'
##
## The following object is masked from 'package:stats':
##
## simulate
grid = setSampling(Terminal = 1, n = 1000)
m1 = setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.var = "x",
time.var = "t", solve.var = "x", xinit = ouFit.ML(wheat)[[2]])
Xwheat = simulate(m1, true.param = list(mu = ouFit.ML(wheat)[[2]], sigma = ouFit.ML(wheat)[[3]],
theta = ouFit.ML(wheat)[[1]]), sampling = grid)
plot(Xwheat)
Increasing the number of simulation to 1000, and time period still is 100 (“yuimu” default: delta = 1/100). Plot the mean of the 1000 silumtion result, different time , very similar result.
simnum = 1000
dist = c(0.31, 0.52, 0.6, 0.7, 0.95)
newsim = function(i) {
simulate(m1, true.param = list(mu = ouFit.ML(wheat)[[2]], sigma = ouFit.ML(wheat)[[3]],
theta = ouFit.ML(wheat)[[1]]))@data@original.data
}
# newsim(1) simulation 1000 times, each time there are 100 time periods
sim = sapply(1:simnum, function(x) newsim(x))
# transfor to time seires format.
m2 = t(sim)
mwheat <- apply(m2, 2, mean)
# plot the mean of the 1000 time simulation for the 100 time periods
plot(mwheat, type = "l")
# find out the quantile to decribe the distribution
tile = sapply(1:100, function(x) quantile(m2[, x], dist))
tile
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## 31% 249.9 248.7 248.3 247.9 247.7 247.4 247.1 246.6 246.3 246.2 246.3
## 52% 249.9 249.9 249.8 250.0 250.1 250.0 250.1 250.4 250.2 250.4 250.1
## 60% 249.9 250.4 250.4 250.7 250.8 251.1 251.2 251.5 251.5 251.6 251.5
## 70% 249.9 251.0 251.3 251.7 252.1 252.4 252.7 253.0 253.0 253.4 253.3
## 95% 249.9 253.4 254.7 255.7 257.1 257.8 258.3 259.1 259.7 260.1 261.1
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 246.0 245.9 245.5 245.3 245.2 245.1 245.1 245.1 244.7 244.8 244.7
## 52% 250.0 250.1 250.1 250.1 250.1 249.9 250.1 250.2 250.2 249.9 250.0
## 60% 251.5 251.6 251.8 251.9 252.1 251.8 251.8 251.9 252.2 252.0 252.2
## 70% 253.6 253.6 253.9 254.2 254.2 254.1 253.9 254.2 254.8 254.7 254.6
## 95% 261.7 262.1 262.0 263.0 263.0 263.9 264.3 264.6 265.5 265.4 266.0
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 244.2 244.2 244.2 243.3 243.3 243.5 242.9 242.8 242.4 242.3 242.2
## 52% 249.6 249.7 249.9 249.6 249.9 249.8 249.6 249.7 249.5 249.6 249.3
## 60% 252.2 252.1 252.2 252.0 252.4 252.4 252.5 251.9 251.9 252.0 251.9
## 70% 254.8 255.4 255.1 255.0 255.4 255.4 255.5 255.7 255.6 255.4 255.3
## 95% 266.4 266.4 266.9 267.3 267.8 268.2 268.5 269.2 270.4 270.9 271.1
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 242.1 242.2 242.1 242.1 242.5 242.4 242.3 242.0 242.0 241.8 241.8
## 52% 249.5 249.5 249.6 249.6 249.7 249.5 249.1 249.3 249.4 249.5 249.5
## 60% 251.9 252.2 252.4 251.7 252.0 252.3 252.1 252.3 252.2 252.0 252.6
## 70% 255.7 255.4 255.5 256.0 255.4 255.9 256.0 256.0 256.1 256.1 256.7
## 95% 271.6 271.6 272.7 272.6 272.4 272.5 272.8 273.6 273.9 273.8 273.5
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 242.2 242.1 241.8 241.6 242.1 241.4 240.9 241.4 241.4 241.1 241.2
## 52% 249.4 249.6 249.8 249.9 250.0 250.0 249.9 249.7 250.2 250.0 250.1
## 60% 252.6 252.7 253.2 253.4 253.3 253.2 253.3 253.4 253.4 253.9 253.5
## 70% 257.1 257.5 257.7 257.7 258.0 257.5 257.2 257.9 257.4 257.9 257.8
## 95% 274.2 274.4 273.7 273.9 273.7 274.1 274.7 274.5 274.7 275.2 276.0
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 241.2 240.5 240.7 240.8 240.3 240.2 240.6 240.4 240.2 240.3 240.2
## 52% 250.1 249.6 249.8 250.1 250.0 249.9 250.7 250.5 250.8 250.2 250.1
## 60% 253.3 253.5 253.3 253.3 253.8 253.6 253.8 254.2 253.6 253.5 253.6
## 70% 258.1 257.6 257.5 257.8 258.5 258.6 258.6 259.0 258.6 258.6 258.6
## 95% 275.9 276.8 276.3 276.2 276.3 276.9 276.7 276.8 277.0 277.7 277.1
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 240.9 240.6 240.6 240.4 240.7 240.4 240.4 240.6 240.5 240.0 240.5
## 52% 250.1 249.7 249.5 249.6 250.1 250.1 250.3 250.6 250.9 250.5 250.5
## 60% 253.5 253.4 253.8 254.1 253.6 254.0 254.0 254.3 254.2 254.2 254.7
## 70% 258.8 258.8 259.0 258.9 259.0 258.7 259.2 259.2 259.5 259.5 259.6
## 95% 278.1 278.3 278.2 279.6 279.4 279.5 279.6 279.7 280.7 280.4 280.8
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 240.8 240.3 240.0 240.5 240.5 240.0 240.3 240.0 240.1 240.4 240.0
## 52% 250.4 250.5 249.9 250.0 250.3 250.7 250.3 250.5 250.8 250.8 250.7
## 60% 254.5 254.5 254.6 254.7 254.5 254.8 254.6 255.4 254.9 255.1 254.6
## 70% 259.7 260.4 260.3 260.3 260.1 260.0 260.5 260.0 260.2 260.8 260.4
## 95% 281.3 280.8 280.7 280.0 281.1 281.4 281.5 281.9 281.7 282.2 282.6
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 240.2 239.9 240.1 239.5 239.8 239.9 240.2 239.3 239.1 239.0 238.6
## 52% 251.0 251.2 251.0 250.9 251.2 251.2 251.6 251.3 250.9 251.2 251.0
## 60% 254.9 255.8 254.9 255.1 255.4 254.9 254.7 255.0 254.6 255.0 254.9
## 70% 260.7 259.7 260.4 260.5 260.6 260.7 260.7 260.9 261.0 260.7 261.2
## 95% 282.5 282.6 282.6 282.5 283.7 284.3 283.3 284.4 283.8 283.4 284.4
## [,100]
## 31% 239.1
## 52% 250.9
## 60% 254.9
## 70% 260.9
## 95% 283.5
One simulation, different time different result. And time period still is 100 (“yuimu” default: delta = 1/100).
require(yuima)
# initial value is 188
m1 = setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.var = "x",
time.var = "t", solve.var = "x", xinit = ouFit.ML(oats)[[2]])
Xoats = simulate(m1, true.param = list(mu = ouFit.ML(oats)[[2]], sigma = ouFit.ML(oats)[[3]],
theta = ouFit.ML(oats)[[1]]))
plot(Xoats)
Increasing the number of simulation to 1000, and time period still is 100 (“yuimu” default: delta = 1/100). Plot the mean of the 1000 silumtion result, different time , very similar result.
simnum = 1000
# specific qunatile (which we can pick any another quantile)
dist = c(0.31, 0.52, 0.6, 0.7, 0.95)
newsim = function(i) {
simulate(m1, true.param = list(mu = ouFit.ML(oats)[[2]], sigma = ouFit.ML(oats)[[3]],
theta = ouFit.ML(oats)[[1]]))@data@original.data
}
# newsim(1) simulation 1000 times, each time there are 100 time periods
sim = sapply(1:simnum, function(x) newsim(x))
# transfor to time seires format.
m2 = t(sim)
moats <- apply(m2, 2, mean)
# plot the mean of the 1000 time simulation for the 100 time periods
plot(moats, type = "l")
# find out the quantile to decribe the distribution
tile = sapply(1:100, function(x) quantile(m2[, x], dist))
tile
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## 31% 186.2 185.8 185.5 185.4 185.2 185.1 185.1 184.9 184.8 184.7 184.6
## 52% 186.2 186.3 186.3 186.2 186.3 186.3 186.3 186.2 186.2 186.3 186.3
## 60% 186.2 186.5 186.6 186.5 186.7 186.7 186.7 186.8 186.9 187.0 186.9
## 70% 186.2 186.7 186.9 187.0 187.2 187.3 187.4 187.4 187.6 187.8 187.7
## 95% 186.2 187.8 188.4 188.9 189.3 189.7 189.9 190.3 190.9 191.1 191.6
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 184.5 184.5 184.4 184.3 184.3 184.3 184.2 184.2 184.0 183.8 183.7
## 52% 186.3 186.3 186.2 186.3 186.3 186.3 186.3 186.3 186.4 186.4 186.5
## 60% 186.9 187.0 186.9 187.0 187.1 187.1 187.1 187.1 187.3 187.2 187.3
## 70% 187.8 187.8 187.9 188.0 188.1 188.2 188.2 188.3 188.3 188.4 188.4
## 95% 191.7 191.8 191.9 192.1 192.4 192.6 193.0 193.1 193.3 193.6 193.7
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 183.7 183.7 183.5 183.6 183.7 183.6 183.5 183.4 183.5 183.4 183.3
## 52% 186.5 186.3 186.2 186.4 186.3 186.3 186.3 186.3 186.5 186.2 186.2
## 60% 187.3 187.2 187.3 187.3 187.3 187.3 187.4 187.5 187.5 187.6 187.4
## 70% 188.5 188.5 188.3 188.5 188.5 188.8 188.7 188.7 188.8 188.8 188.8
## 95% 193.9 194.1 194.8 194.7 194.5 194.9 195.0 195.0 194.9 195.1 195.6
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 183.3 183.1 183.1 183.2 183.1 182.7 182.6 182.5 182.4 182.6 182.6
## 52% 186.3 186.2 186.2 186.2 186.2 186.1 186.1 186.1 186.1 186.2 186.3
## 60% 187.5 187.3 187.3 187.4 187.3 187.4 187.4 187.5 187.5 187.5 187.4
## 70% 188.8 188.8 188.9 188.9 189.2 189.2 189.2 189.2 189.1 189.1 189.1
## 95% 195.5 195.4 195.5 195.5 195.6 195.7 195.9 196.2 196.0 196.0 196.3
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 182.8 182.7 182.6 182.7 182.6 182.6 182.8 182.6 182.5 182.4 182.5
## 52% 186.1 186.1 186.2 186.0 186.1 186.2 186.3 186.3 186.3 186.4 186.3
## 60% 187.4 187.2 187.3 187.3 187.4 187.4 187.6 187.7 187.7 187.7 187.8
## 70% 189.1 189.3 189.4 189.3 189.3 189.4 189.3 189.4 189.5 189.7 189.7
## 95% 196.1 196.4 196.2 196.7 196.7 196.8 196.8 196.7 197.1 197.2 197.5
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 182.6 182.3 182.3 182.3 182.4 182.2 182.1 181.9 181.9 181.9 181.9
## 52% 186.4 186.3 186.5 186.3 186.4 186.3 186.2 185.9 186.2 186.3 186.3
## 60% 187.7 187.9 187.8 187.8 187.8 187.8 187.8 187.7 187.8 187.8 187.9
## 70% 189.6 189.6 189.7 189.8 189.9 189.9 189.7 190.0 189.9 190.0 190.0
## 95% 197.6 197.8 198.1 198.1 198.3 198.2 198.3 198.3 198.4 198.6 198.5
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 181.9 181.9 182.0 181.9 181.8 182.0 182.0 182.1 182.1 182.2 182.1
## 52% 186.2 186.3 186.3 186.5 186.6 186.6 186.5 186.6 186.4 186.5 186.5
## 60% 187.8 187.6 187.9 188.1 188.1 188.1 188.1 188.2 188.1 188.1 188.1
## 70% 189.7 189.8 189.9 190.0 190.0 190.2 190.0 190.0 190.0 190.1 190.3
## 95% 198.7 198.8 199.0 199.0 199.0 199.3 199.5 200.1 199.8 199.7 199.9
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 182.2 182.3 182.2 181.9 181.7 181.7 181.4 181.3 181.4 181.5 181.4
## 52% 186.4 186.6 186.6 186.5 186.6 186.7 186.8 186.6 186.6 186.7 186.7
## 60% 188.1 188.3 188.1 188.2 188.1 188.4 188.3 188.2 188.3 188.2 188.3
## 70% 190.1 190.1 190.3 190.4 190.4 190.6 190.4 190.3 190.5 190.6 190.5
## 95% 199.9 200.4 199.9 200.3 200.6 200.2 200.2 200.0 199.9 200.0 200.0
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 181.3 181.4 181.3 181.2 180.9 180.9 181.2 181.2 181.2 181.1 181.2
## 52% 186.6 186.9 186.8 186.7 186.6 186.6 186.5 186.4 186.5 186.3 186.5
## 60% 188.4 188.3 188.4 188.4 188.5 188.5 188.4 188.3 188.4 188.4 188.3
## 70% 190.5 190.6 190.7 190.8 190.7 190.7 190.7 190.7 190.8 190.9 190.9
## 95% 200.0 200.2 200.5 200.6 200.7 200.8 201.2 201.0 200.9 201.0 200.8
## [,100]
## 31% 181.1
## 52% 186.5
## 60% 188.3
## 70% 191.0
## 95% 200.9
One simulation, different time different result. And time period still is 100 (“yuimu” default: delta = 1/100).
require(yuima)
# initial value is 188
m1 = setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.var = "x",
time.var = "t", solve.var = "x", xinit = ouFit.ML(barley)[[2]])
Xbarley = simulate(m1, true.param = list(mu = ouFit.ML(barley)[[2]], sigma = ouFit.ML(barley)[[3]],
theta = ouFit.ML(barley)[[1]]))
plot(Xbarley)
Increasing the number of simulation to 1000, and time period still is 100 (“yuimu” default: delta = 1/100). Plot the mean of the 1000 silumtion result, different time , very similar result.
simnum = 1000
# specific qunatile (which we can pick any another quantile)
dist = c(0.31, 0.52, 0.6, 0.7, 0.95)
newsim = function(i) {
simulate(m1, true.param = list(mu = ouFit.ML(barley)[[2]], sigma = ouFit.ML(barley)[[3]],
theta = ouFit.ML(barley)[[1]]))@data@original.data
}
# newsim(1) simulation 1000 times, each time there are 100 time periods
sim = sapply(1:simnum, function(x) newsim(x))
# transfor to time seires format.
m2 = t(sim)
mbarley <- apply(m2, 2, mean)
# plot the mean of the 1000 time simulation for the 100 time periods
plot(mbarley, type = "l")
# find out the quantile to decribe the distribution
tile = sapply(1:100, function(x) quantile(m2[, x], dist))
tile
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## 31% 188.8 188.2 188.0 187.7 187.5 187.4 187.3 187.2 187.1 187.0 187.0
## 52% 188.8 188.8 188.8 188.8 188.9 188.8 189.0 188.9 188.8 188.8 188.9
## 60% 188.8 189.1 189.1 189.2 189.3 189.3 189.5 189.4 189.5 189.5 189.5
## 70% 188.8 189.4 189.6 189.7 189.9 190.0 190.2 190.3 190.3 190.4 190.5
## 95% 188.8 190.7 191.3 191.8 192.4 192.7 193.2 193.6 194.0 194.3 194.4
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 187.0 187.1 187.0 186.9 186.9 186.7 186.8 186.5 186.5 186.3 186.2
## 52% 188.9 189.1 189.1 189.1 189.0 188.9 188.9 189.1 189.0 189.1 189.0
## 60% 189.6 189.7 189.7 189.9 189.9 189.9 189.9 190.0 190.1 190.1 190.1
## 70% 190.6 190.7 190.8 190.9 191.0 191.1 191.3 191.3 191.4 191.4 191.6
## 95% 194.7 195.0 195.4 195.6 195.9 195.9 196.4 196.6 196.7 197.3 197.4
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 186.2 186.2 186.0 186.0 186.0 186.0 185.8 185.9 185.8 185.8 185.6
## 52% 189.1 189.2 189.0 188.9 189.2 189.2 189.3 189.2 189.1 189.0 189.1
## 60% 190.3 190.2 190.2 190.2 190.2 190.2 190.1 190.4 190.4 190.3 190.4
## 70% 191.5 191.6 191.7 191.8 191.8 191.7 191.7 191.8 192.0 191.8 191.9
## 95% 197.5 197.3 197.8 198.0 198.1 198.3 198.2 198.7 198.9 199.0 198.9
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 185.5 185.8 185.7 185.4 185.5 185.1 185.0 185.0 185.0 185.0 185.1
## 52% 189.1 189.2 189.3 189.2 189.2 189.1 189.1 189.2 189.3 189.0 189.1
## 60% 190.4 190.4 190.4 190.4 190.4 190.3 190.5 190.2 190.5 190.3 190.4
## 70% 192.2 192.1 192.1 192.2 192.0 192.1 191.9 192.0 192.1 192.3 192.4
## 95% 199.1 199.1 199.2 199.6 199.4 199.7 199.4 199.6 199.8 199.9 200.4
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 185.0 184.7 184.9 184.9 184.7 184.6 184.8 184.6 184.6 184.5 184.4
## 52% 188.9 189.1 189.2 189.0 188.8 188.9 188.9 188.9 188.7 188.9 189.1
## 60% 190.5 190.6 190.5 190.5 190.4 190.4 190.4 190.4 190.5 190.7 190.8
## 70% 192.7 192.7 192.5 192.4 192.5 192.5 192.4 192.6 192.5 192.9 192.8
## 95% 200.6 200.8 200.7 201.0 201.1 201.7 201.4 201.6 201.5 201.7 201.5
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 184.4 184.2 184.2 184.2 184.0 184.2 183.9 183.9 183.7 183.9 183.8
## 52% 188.8 188.9 188.7 188.8 188.5 188.8 188.7 188.4 188.5 188.6 188.5
## 60% 190.8 190.9 190.7 190.8 190.8 191.0 190.6 190.5 190.4 190.3 190.5
## 70% 192.9 193.2 193.1 192.9 193.0 192.9 193.2 193.0 193.1 193.0 193.0
## 95% 201.9 202.3 202.3 202.5 202.4 202.7 202.8 203.0 203.4 203.2 203.2
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 183.7 183.8 183.9 183.8 183.9 183.9 184.0 183.8 183.9 183.7 183.8
## 52% 188.5 188.4 188.6 188.6 188.6 188.5 188.9 189.0 188.8 188.6 188.5
## 60% 190.5 190.6 190.7 190.6 190.5 190.5 190.5 190.6 190.6 190.5 190.6
## 70% 193.1 193.0 192.9 192.7 193.2 192.9 192.9 192.9 193.0 193.1 193.3
## 95% 203.9 203.8 204.2 203.8 203.7 203.8 204.2 204.3 204.0 203.9 204.1
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 183.8 183.7 183.4 183.5 183.5 183.7 183.7 183.7 183.7 183.7 183.5
## 52% 188.7 188.9 188.9 189.0 188.9 188.9 188.9 188.8 188.5 188.9 188.6
## 60% 190.6 190.5 190.7 190.8 190.7 190.7 190.6 190.7 190.6 191.0 191.1
## 70% 193.3 193.3 193.4 193.6 193.5 193.6 193.5 193.5 193.8 194.0 193.9
## 95% 204.1 204.0 203.8 204.4 204.6 205.2 205.1 205.2 205.3 205.0 205.2
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 183.6 183.4 183.2 183.1 183.0 183.0 183.1 183.1 183.3 183.4 183.2
## 52% 188.8 188.9 188.9 188.9 188.9 188.8 189.0 189.1 189.2 189.5 189.5
## 60% 191.2 191.1 191.0 190.9 190.9 191.2 191.0 191.1 191.5 191.4 191.4
## 70% 194.0 194.0 193.9 193.9 193.9 194.1 194.2 194.4 194.3 194.4 194.5
## 95% 205.4 205.5 205.6 206.0 206.1 206.0 206.1 205.9 205.9 206.2 206.4
## [,100]
## 31% 183.1
## 52% 189.4
## 60% 191.5
## 70% 194.8
## 95% 206.6
One simulation, different time different result. And time period still is 100 (“yuimu” default: delta = 1/100).
require(yuima)
# initial value is 188
m1 = setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.var = "x",
time.var = "t", solve.var = "x", xinit = ouFit.ML(canola)[[2]])
Xcanola = simulate(m1, true.param = list(mu = ouFit.ML(canola)[[2]], sigma = ouFit.ML(canola)[[3]],
theta = ouFit.ML(canola)[[1]]))
plot(Xcanola)
Increasing the number of simulation to 1000, and time period still is 100 (“yuimu” default: delta = 1/100). Plot the mean of the 1000 silumtion result, different time , very similar result.
simnum = 1000
# specific qunatile (which we can pick any another quantile)
dist = c(0.31, 0.52, 0.6, 0.7, 0.95)
newsim = function(i) {
simulate(m1, true.param = list(mu = ouFit.ML(canola)[[2]], sigma = ouFit.ML(canola)[[3]],
theta = ouFit.ML(canola)[[1]]))@data@original.data
}
# newsim(1) simulation 1000 times, each time there are 100 time periods
sim = sapply(1:simnum, function(x) newsim(x))
# transfor to time seires format.
m2 = t(sim)
mcanola <- apply(m2, 2, mean)
# plot the mean of the 1000 time simulation for the 100 time periods
plot(mcanola, type = "l")
# find out the quantile to decribe the distribution
tile = sapply(1:100, function(x) quantile(m2[, x], dist))
tile
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## 31% 507.7 506.7 506.3 505.8 505.7 505.6 505.2 505.0 505.0 504.7 504.4
## 52% 507.7 507.8 507.6 507.7 507.6 507.6 507.6 507.8 507.8 507.9 507.9
## 60% 507.7 508.1 508.2 508.3 508.4 508.4 508.4 508.7 508.7 508.8 508.9
## 70% 507.7 508.6 508.8 509.1 509.4 509.5 509.7 510.0 510.1 510.3 510.2
## 95% 507.7 510.5 511.9 512.6 513.3 514.4 514.6 515.2 515.9 516.3 516.6
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 504.3 504.3 504.2 504.1 503.6 503.8 503.7 503.6 503.1 502.9 502.8
## 52% 507.4 507.5 507.5 507.4 507.5 507.6 507.6 507.6 507.4 507.6 507.1
## 60% 508.8 508.6 508.6 508.8 508.8 508.9 509.0 509.1 509.0 509.2 509.1
## 70% 510.2 510.4 510.4 510.5 510.6 510.7 510.8 510.7 511.3 511.3 511.3
## 95% 517.2 517.4 517.9 518.1 518.3 518.1 518.5 518.6 519.0 519.1 519.8
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 502.7 502.6 502.2 502.2 502.3 502.0 501.9 501.7 501.9 501.5 501.4
## 52% 507.4 507.2 507.4 507.5 507.4 507.3 507.3 507.3 507.2 507.3 507.4
## 60% 509.1 509.1 509.1 509.0 508.9 509.3 509.3 509.4 509.1 509.3 509.2
## 70% 511.4 511.2 510.9 511.1 511.2 511.3 511.6 511.7 511.8 511.8 512.0
## 95% 519.8 520.6 521.1 521.3 521.6 521.5 522.0 522.6 522.7 522.7 523.3
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 501.5 501.5 501.4 501.5 501.3 501.5 501.3 501.2 501.3 501.1 501.3
## 52% 507.8 507.5 507.4 507.7 507.4 507.1 507.2 507.5 507.4 507.3 507.2
## 60% 509.3 509.3 509.5 509.8 509.6 509.5 509.5 509.6 509.4 509.5 509.7
## 70% 512.2 512.2 512.1 512.1 512.2 512.4 512.5 512.6 513.1 512.9 512.9
## 95% 523.6 524.4 524.1 524.2 524.5 524.7 524.8 524.9 525.1 525.4 526.0
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 500.9 501.0 500.8 500.5 500.8 500.9 500.9 500.6 500.5 500.6 500.5
## 52% 507.2 507.6 507.7 507.6 507.5 507.4 507.4 507.5 507.5 507.5 507.6
## 60% 510.0 510.3 510.2 510.3 510.2 510.2 510.0 509.9 510.1 510.0 510.0
## 70% 513.0 513.4 513.4 513.4 513.6 513.6 513.6 513.5 513.9 513.8 513.9
## 95% 526.2 526.2 526.3 526.5 526.8 526.8 526.9 527.6 527.8 527.5 528.2
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 500.5 500.8 500.8 500.6 500.7 500.6 500.4 500.4 500.2 499.7 500.2
## 52% 507.6 507.6 507.8 507.7 507.6 507.8 507.7 507.6 507.6 507.3 507.6
## 60% 510.1 510.0 510.2 510.4 510.4 510.4 510.5 510.5 510.6 510.5 510.6
## 70% 513.7 513.9 514.4 514.3 514.4 514.1 513.9 513.8 513.9 514.2 514.3
## 95% 528.7 528.7 528.4 528.2 528.4 529.6 529.2 529.9 530.2 529.6 530.1
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 500.1 500.0 500.1 500.1 500.1 499.8 499.3 499.7 499.5 499.4 499.0
## 52% 507.6 507.3 507.5 507.5 507.0 507.1 507.4 507.3 507.5 507.5 507.3
## 60% 510.7 510.3 510.4 510.5 510.5 510.4 510.6 510.3 510.1 510.5 510.7
## 70% 514.4 514.5 514.7 514.8 514.8 514.9 514.9 514.7 514.7 514.7 514.9
## 95% 531.0 530.6 531.4 531.1 531.3 531.6 532.3 532.6 532.3 532.4 532.7
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 498.9 498.8 498.8 499.0 498.8 499.0 499.0 498.8 499.0 499.2 498.8
## 52% 507.6 507.6 507.7 507.4 507.4 507.9 507.7 508.0 507.5 507.7 507.3
## 60% 510.8 510.9 510.6 510.4 510.4 510.4 510.3 510.4 510.5 510.4 510.5
## 70% 515.3 515.2 514.7 514.5 514.7 514.7 514.9 514.6 514.9 514.4 514.4
## 95% 533.3 533.2 533.5 533.5 532.6 532.4 533.0 533.6 533.4 533.4 533.2
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 498.8 498.7 498.9 498.8 499.1 498.9 499.4 499.3 498.9 498.4 498.3
## 52% 507.3 507.3 507.5 507.4 507.7 507.4 507.2 507.4 507.0 506.9 507.2
## 60% 510.6 510.5 510.5 510.8 510.6 510.7 510.7 510.5 510.3 510.1 510.1
## 70% 514.7 514.8 514.6 515.1 514.9 515.0 515.0 515.1 515.8 515.5 515.7
## 95% 533.9 533.0 533.5 534.3 534.1 534.9 535.1 534.8 535.7 535.4 535.5
## [,100]
## 31% 498.4
## 52% 507.6
## 60% 510.5
## 70% 516.0
## 95% 534.9
One simulation, different time different result. And time period still is 100 (“yuimu” default: delta = 1/100).
require(yuima)
# initial value is 188
m1 = setModel(drift = "theta*(mu-x)", diffusion = "sigma", state.var = "x",
time.var = "t", solve.var = "x", xinit = ouFit.ML(dry.peas)[[2]])
Xdry.peas = simulate(m1, true.param = list(mu = ouFit.ML(dry.peas)[[2]], sigma = ouFit.ML(dry.peas)[[3]],
theta = ouFit.ML(dry.peas)[[1]]))
plot(Xdry.peas)
Increasing the number of simulation to 1000, and time period still is 100 (“yuimu” default: delta = 1/100). Plot the mean of the 1000 silumtion result, different time , very similar result.
simnum = 1000
# specific qunatile (which we can pick any another quantile)
dist = c(0.31, 0.52, 0.6, 0.7, 0.95)
newsim = function(i) {
simulate(m1, true.param = list(mu = ouFit.ML(dry.peas)[[2]], sigma = ouFit.ML(dry.peas)[[3]],
theta = ouFit.ML(dry.peas)[[1]]))@data@original.data
}
# newsim(1) simulation 1000 times, each time there are 100 time periods
sim = sapply(1:simnum, function(x) newsim(x))
# transfor to time seires format.
m2 = t(sim)
mdry.peas <- apply(m2, 2, mean)
# plot the mean of the 1000 time simulation for the 100 time periods
plot(mdry.peas, type = "l")
# find out the quantile to decribe the distribution
tile = sapply(1:100, function(x) quantile(m2[, x], dist))
tile
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
## 31% 279.5 278.8 278.6 278.4 278.1 277.9 277.8 277.7 277.4 277.3 277.6
## 52% 279.5 279.5 279.6 279.6 279.6 279.6 279.7 279.7 279.6 279.8 279.7
## 60% 279.5 279.8 280.1 280.2 280.3 280.3 280.4 280.4 280.6 280.6 280.5
## 70% 279.5 280.2 280.5 280.8 281.0 281.1 281.3 281.5 281.7 281.6 281.8
## 95% 279.5 281.9 282.9 283.4 284.1 284.9 285.5 286.1 286.2 286.2 286.9
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 277.3 277.0 276.8 276.9 276.9 276.9 276.7 276.6 276.4 276.4 276.2
## 52% 279.6 279.7 279.6 279.8 279.9 279.8 279.9 279.9 279.9 279.8 280.0
## 60% 280.5 280.6 280.8 281.0 281.1 281.0 280.9 281.0 281.2 281.3 281.5
## 70% 281.8 282.0 282.0 282.3 282.5 282.5 282.7 282.7 282.8 283.0 283.2
## 95% 287.1 287.5 287.6 288.0 287.9 288.6 288.4 288.4 288.8 289.8 289.8
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 276.3 276.2 276.2 276.2 276.0 276.1 275.9 275.7 275.6 275.5 275.8
## 52% 280.1 280.1 280.2 279.7 279.7 279.8 279.9 280.1 280.0 280.3 280.2
## 60% 281.4 281.5 281.7 281.6 281.5 281.5 281.7 281.7 281.8 281.9 282.0
## 70% 283.4 283.5 283.6 283.6 283.8 284.0 283.7 283.9 284.1 284.3 284.2
## 95% 290.2 290.8 291.3 291.7 292.2 292.0 292.0 292.1 292.4 292.7 293.1
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 275.3 275.5 275.3 275.3 275.2 275.3 275.4 275.4 275.2 275.5 274.9
## 52% 280.1 280.5 280.4 280.3 280.4 280.6 280.8 280.5 280.5 280.4 280.2
## 60% 282.2 282.1 282.0 281.9 282.0 282.1 282.1 282.3 282.2 282.2 282.2
## 70% 284.2 284.3 284.5 284.6 284.4 284.5 284.5 284.9 284.6 284.5 284.6
## 95% 292.9 292.9 293.2 293.4 293.9 294.4 294.7 294.6 294.9 295.3 295.6
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 275.1 275.0 274.8 274.8 274.5 274.3 274.5 274.6 274.5 274.3 274.7
## 52% 280.3 280.2 280.4 280.3 280.3 280.3 280.2 280.1 279.8 280.2 280.0
## 60% 282.5 282.3 282.4 282.3 282.2 282.3 282.4 282.5 282.4 282.4 282.6
## 70% 285.0 284.9 284.8 284.7 284.7 285.0 285.1 285.4 285.4 285.3 285.3
## 95% 295.5 295.5 295.5 295.8 296.6 296.0 296.6 296.5 296.7 297.0 297.4
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 274.4 274.3 274.1 274.3 274.1 274.3 274.4 274.5 274.0 274.0 274.3
## 52% 280.0 280.1 280.0 280.1 280.3 280.2 280.5 280.1 280.4 280.3 280.2
## 60% 282.5 282.3 282.4 282.4 282.5 282.7 282.4 282.5 282.7 282.7 282.9
## 70% 285.4 285.5 285.3 285.5 285.7 285.7 285.8 286.1 286.0 285.9 286.1
## 95% 297.6 297.8 298.2 298.0 298.1 297.9 298.5 297.7 298.6 298.3 298.8
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 273.9 274.0 273.7 273.5 273.5 273.8 273.9 273.9 273.7 273.8 273.7
## 52% 280.5 280.3 280.0 280.2 280.2 280.4 280.7 280.7 280.5 280.8 280.7
## 60% 282.6 282.7 282.8 282.7 282.9 282.7 282.9 283.1 283.1 283.2 283.4
## 70% 286.2 286.2 286.2 286.1 286.7 286.7 286.8 286.4 286.4 286.5 286.5
## 95% 298.7 298.7 298.8 298.5 298.9 299.3 299.7 299.5 299.8 300.0 300.4
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 273.7 273.8 273.7 273.3 273.3 273.1 273.3 273.3 273.4 273.6 273.3
## 52% 280.4 280.2 280.3 280.5 280.1 280.3 280.4 280.5 280.5 280.2 280.2
## 60% 283.3 283.4 283.2 283.3 283.1 283.1 283.1 282.8 282.8 283.0 283.2
## 70% 286.7 286.8 287.0 286.9 286.6 286.6 287.0 286.9 287.1 286.7 287.0
## 95% 301.0 300.6 301.3 301.1 301.3 300.9 301.0 300.7 300.9 300.8 301.2
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 273.5 273.7 273.6 273.8 273.5 273.2 273.0 272.9 273.0 273.1 272.9
## 52% 280.6 280.4 280.6 280.3 280.3 280.7 280.6 280.8 280.4 280.3 280.2
## 60% 283.6 283.4 283.5 283.4 283.3 283.2 283.6 283.4 283.4 283.5 283.5
## 70% 287.3 287.6 287.5 287.5 288.0 287.7 287.5 287.1 287.6 287.5 287.4
## 95% 301.4 301.0 301.3 301.7 301.4 301.5 301.5 301.9 302.0 302.5 302.4
## [,100]
## 31% 272.7
## 52% 280.5
## 60% 283.7
## 70% 287.5
## 95% 302.5
Recall “sde.sim” package uses this notation:
\[ dX_t = (\theta_1 - \theta_2X_t)dt +\theta_3dW_t, X_0 = x_0 \]
For “sde.sim” packages, each price simulation uses different random seed. When we try the same random seed, all prices have the same trend.
set.seed(123)
thetaw <- c(ouFit.ML(wheat)[[1]] * ouFit.ML(wheat)[[2]], ouFit.ML(wheat)[[1]],
ouFit.ML(wheat)[[3]])
Xwheat <- sde.sim(X0 = ouFit.ML(wheat)[[2]], model = "OU", theta = thetaw, N = 1000,
delta = 1)
## Error: could not find function "sde.sim"
plot(Xwheat, main = " Ornstein - Uhlenbeck ")
set.seed(321)
thetao <- c(ouFit.ML(oats)[[1]] * ouFit.ML(oats)[[2]], ouFit.ML(oats)[[1]],
ouFit.ML(oats)[[3]])
Xoats <- sde.sim(X0 = ouFit.ML(oats)[[2]], model = "OU", theta = thetao, N = 1000,
delta = 1)
## Error: could not find function "sde.sim"
plot(Xoats, main = " Ornstein - Uhlenbeck ")
set.seed(132)
thetab <- c(ouFit.ML(barley)[[1]] * ouFit.ML(barley)[[2]], ouFit.ML(barley)[[1]],
ouFit.ML(barley)[[3]])
# set initial value as the long term mean mu.
Xbarley <- sde.sim(X0 = ouFit.ML(barley)[[2]], model = "OU", theta = thetab,
N = 1000, delta = 1)
## Error: could not find function "sde.sim"
plot(Xbarley, main = " Ornstein - Uhlenbeck ")
set.seed(312)
thetac <- c(ouFit.ML(canola)[[1]] * ouFit.ML(canola)[[2]], ouFit.ML(canola)[[1]],
ouFit.ML(canola)[[3]])
Xcanola <- sde.sim(X0 = ouFit.ML(canola)[[2]], model = "OU", theta = thetac,
N = 1000, delta = 1)
## Error: could not find function "sde.sim"
plot(Xcanola, main = " Ornstein - Uhlenbeck ")
set.seed(231)
thetad <- c(ouFit.ML(dry.peas)[[1]] * ouFit.ML(dry.peas)[[2]], ouFit.ML(dry.peas)[[1]],
ouFit.ML(dry.peas)[[3]])
Xdry.peas <- sde.sim(X0 = ouFit.ML(dry.peas)[[2]], model = "OU", theta = thetad,
N = 1000, delta = 1)
## Error: could not find function "sde.sim"
plot(Xdry.peas, main = " Ornstein - Uhlenbeck ")
This documents show we can use ouFit.ML function to estimate the parameter \( \theta \), \( \mu \), and \( \sigma \). And them use the “yuima” package to generate the simulative price vector/matrices.