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 a stochastic process which is stationary, Gaussian, and Markovian.
Over time, the process tends to drift towards its long-term mean: such a process is called mean-reverting.
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.
“The Ornstein-Uhlenbeck process is one of several approaches used to model (with modifications) interest rates, currency exchange rates, and commodity prices stochastically. The parameter \( \mu \) represents the equilibrium or mean value supported by fundamentals; \( \sigma \) the degree of volatility around it caused by shocks, and \( \theta \) the rate by which these shocks dissipate and the variable reverts towards the mean.”(Wikipedia)
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.
In excel, we transform the date format to month/day/year, which save some time.
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.
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') # not necessary
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
paraWheat <- c(ouFit.ML(wheat)[[1]], ouFit.ML(wheat)[[2]], ouFit.ML(wheat)[[3]])
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
paraOats <- c(ouFit.ML(oats)[[1]], ouFit.ML(oats)[[2]], ouFit.ML(oats)[[3]])
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
paraBarley <- c(ouFit.ML(barley)[[1]], ouFit.ML(barley)[[2]], ouFit.ML(barley)[[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
paraCanola <- c(ouFit.ML(canola)[[1]], ouFit.ML(canola)[[2]], ouFit.ML(canola)[[3]])
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
paraDrypeas <- c(ouFit.ML(dry.peas)[[1]], ouFit.ML(dry.peas)[[2]], ouFit.ML(dry.peas)[[3]])
para <- data.frame(rbind(paraWheat, paraOats, paraBarley, paraCanola, paraDrypeas))
rownames(para) <- c("Wheat", "Oats", "Barley", "Canola", "Drypeas")
colnames(para) <- c("Theta", "Mu", "Delta")
para
## Theta Mu Delta
## Wheat 0.11982 249.9 21.993
## Oats 0.08256 186.2 9.471
## Barley 0.04971 188.8 11.290
## Canola 0.04689 507.7 18.465
## Drypeas 0.05435 279.5 14.076
plot(wheat, type = "l", lwd = 2.5, ylim = c(120, 650), xlab = "Crop prices",
col = 1)
lines(oats, type = "l", lwd = 2.5, col = 2)
lines(barley, type = "l", lwd = 2.5, col = 3)
lines(canola, type = "l", lwd = 2.5, col = 4)
lines(dry.peas, type = "l", lwd = 2.5, col = 5)
legend("topleft", c("wheat", "oats", "barley", "canola", "dry.peas"), lty = c(1,
1), lwd = c(2.5, 2.5), col = c(1, 2, 3, 4, 5))
There are two packages we can use to generate Vasicek process.
The following simulation equation can be used for generating paths (sampled with fixed time steps of \( \delta \)). The equation is an exact solution of the SDE. \[ S_{i+1} = S+i e^{-\theta \delta} + \mu (1- e^{-\theta \delta}) + \delta \sqrt(\frac{1- e^{-2 \theta \delta}}{2 \theta}) N_{0,1} \]
# Brownian motion
# Example from Dixit and Pindyck, 1994, pp.74-76
# Simple mean-reverting process:% dx = nu (xbar - x) dt + sigma dz
#
OU.sim <- function(para = c(0.01, 200, 1), x0= 1, periods=100){
periods=100; #Number of periods
nu = c(0.0, 0.01, 0.02, 0.5); # speed of reversion
sigma = ouFit.ML(wheat)[[2]]; # in monthly terms
ones = c(1,1,1,1) # all one vector
sigma2 = (((sigma^2)*(2*nu)).*(ones-exp(-2*nu))).*0.5; #dt=1;
xbar=c(1, 1, 1, 1); # Level to which x tends to revert, or normal level
x=matrix(rep(0, periods*4), periods, 4) # all zero, 100*4 matrix
epsilon=matrix(rnorm(periods*4),periods,4); # random normal number, 100*4 matrix
x[1,]=c(1,1,1,1); # Starting value of first row of x
i=2;
while i<periods+1
# x(i, :)=x(i-1, :)+nu.*(xbar-x(i-1, :)) + sigma*epsilon(i-1, :);
x(i, :) = x(i-1, :) + xbar.*(ones(1,4)-exp(-nu)) + x(i-1, :).* ...
(exp(-nu)-ones(1,4)) + sigma*epsilon(i-1, :);
i=i+1;
end;
#figure (1)
#plot(x);
#xlabel('Period');
#ylabel('Value of x');
#title('Ornstein-Uhlenbeck Stochastic Mean-Reversion Process');
#legend('0.00', '0.01', '0.02','0.50');
}
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 delta of t
\[ dt=(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.9 248.4 247.9 247.9 247.6 247.4 247.1 246.9 247.0 246.8
## 52% 249.9 250.1 250.1 250.0 250.2 250.4 250.6 250.4 250.3 250.4 250.5
## 60% 249.9 250.5 250.9 251.0 251.2 251.5 251.7 251.8 251.8 251.9 251.9
## 70% 249.9 251.1 251.7 251.9 252.5 252.7 253.1 253.2 253.6 253.6 253.6
## 95% 249.9 253.5 255.3 256.5 257.2 258.0 258.6 259.1 259.9 260.4 260.6
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 246.9 246.7 246.4 246.4 245.9 245.8 245.6 245.6 245.4 245.1 245.2
## 52% 250.4 250.4 250.6 250.6 250.5 250.5 250.6 250.7 250.7 250.4 250.6
## 60% 251.9 251.7 252.1 252.0 252.1 252.1 252.4 252.6 252.4 252.2 252.5
## 70% 253.7 253.6 254.1 254.2 253.9 254.1 254.5 254.8 255.0 255.0 254.9
## 95% 261.6 262.1 262.5 262.8 263.1 263.4 264.2 264.4 264.7 264.7 265.0
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 244.9 244.7 244.8 244.9 244.4 244.4 244.5 244.2 244.0 243.6 243.5
## 52% 250.4 250.3 250.5 250.1 250.3 250.3 250.2 250.1 250.4 250.3 250.6
## 60% 252.3 252.3 252.4 252.4 252.4 252.3 252.7 252.4 252.7 252.8 252.9
## 70% 255.3 255.4 255.5 255.7 255.6 255.7 256.1 255.8 255.7 255.8 255.7
## 95% 265.4 266.2 266.9 266.8 267.0 267.5 268.1 268.4 268.5 269.3 269.3
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 243.2 243.3 243.2 242.9 242.9 242.7 243.1 242.9 242.5 242.3 242.1
## 52% 250.6 250.8 250.7 250.7 250.5 250.5 250.2 250.2 250.4 249.7 250.5
## 60% 252.9 252.9 252.9 253.0 252.9 253.2 253.0 252.9 252.9 253.0 252.9
## 70% 256.1 256.1 256.3 256.4 256.5 256.8 256.9 256.9 256.7 256.2 256.7
## 95% 269.1 269.7 270.2 270.3 270.4 269.7 269.6 270.9 271.6 272.3 272.1
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 241.7 241.9 242.1 241.9 242.3 241.9 242.0 242.2 241.8 241.8 241.5
## 52% 250.0 250.4 250.5 250.0 250.2 250.2 249.9 250.0 250.0 250.3 250.3
## 60% 253.4 253.2 253.0 253.1 253.4 253.7 253.3 252.9 252.9 253.0 253.1
## 70% 257.0 256.7 257.2 257.2 257.2 257.3 257.7 257.5 257.1 256.8 257.5
## 95% 272.0 272.5 272.8 272.9 273.1 273.6 273.9 273.0 272.6 273.5 273.6
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 241.4 241.2 241.3 240.8 240.9 240.7 241.1 241.3 241.4 240.9 241.4
## 52% 250.9 250.9 250.6 250.9 251.1 250.5 250.9 250.9 250.9 251.0 250.7
## 60% 253.6 253.8 254.0 254.2 253.8 254.2 253.9 254.1 254.1 254.0 254.2
## 70% 257.2 257.4 257.6 257.2 258.1 258.3 258.7 259.2 258.8 259.0 258.7
## 95% 274.5 273.9 275.4 276.2 275.7 275.8 275.9 276.6 277.0 277.1 276.5
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 241.3 241.3 241.2 241.2 241.0 241.0 240.8 240.8 240.5 240.4 240.5
## 52% 250.6 250.9 250.3 250.3 250.0 250.3 250.6 250.6 250.3 250.3 250.0
## 60% 254.3 254.1 253.9 253.8 253.7 254.2 253.9 254.0 253.8 253.8 253.7
## 70% 258.8 258.8 258.7 258.0 258.1 257.9 258.1 258.2 258.5 258.1 258.2
## 95% 276.9 277.1 276.7 276.2 277.8 277.7 277.7 278.3 278.7 278.8 278.5
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 240.4 240.8 240.9 240.4 240.7 240.2 240.3 239.9 239.6 239.3 239.4
## 52% 250.1 250.2 250.3 250.4 250.0 250.0 250.1 250.0 249.4 249.8 250.0
## 60% 254.0 253.6 253.4 253.6 253.9 254.1 254.2 253.7 253.9 253.8 253.8
## 70% 258.7 258.9 259.2 258.7 258.4 258.7 258.8 258.9 259.2 259.0 258.7
## 95% 279.7 280.2 280.3 280.2 280.5 281.6 282.0 281.9 281.6 281.9 282.4
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 239.7 239.9 239.7 239.8 240.1 239.9 239.3 239.1 239.1 239.2 239.3
## 52% 249.6 249.9 249.9 249.8 249.9 249.8 250.3 250.4 250.3 250.3 250.2
## 60% 253.3 253.0 253.7 253.4 254.1 253.8 253.8 254.2 254.0 253.9 254.3
## 70% 259.1 259.4 259.2 259.1 259.4 259.4 259.4 259.5 259.5 259.7 259.6
## 95% 282.2 282.0 282.3 282.7 282.5 282.4 282.7 282.5 282.0 283.0 282.4
## [,100]
## 31% 239.3
## 52% 249.9
## 60% 254.1
## 70% 260.1
## 95% 282.0
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.6 185.5 185.4 185.4 185.2 185.1 185.0 184.8 184.7
## 52% 186.2 186.4 186.3 186.4 186.4 186.4 186.4 186.4 186.3 186.4 186.4
## 60% 186.2 186.5 186.6 186.7 186.8 186.9 186.9 186.9 186.9 186.9 187.0
## 70% 186.2 186.8 187.0 187.2 187.3 187.5 187.5 187.5 187.5 187.6 187.7
## 95% 186.2 187.8 188.6 189.0 189.3 189.7 190.1 190.4 190.6 191.0 191.0
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 184.7 184.6 184.6 184.5 184.4 184.4 184.4 184.3 184.2 184.1 183.9
## 52% 186.4 186.4 186.4 186.5 186.5 186.5 186.3 186.6 186.6 186.5 186.5
## 60% 187.1 187.0 187.1 187.2 187.2 187.2 187.2 187.4 187.4 187.4 187.3
## 70% 187.8 187.9 187.9 188.0 188.2 188.2 188.4 188.3 188.4 188.4 188.4
## 95% 191.2 191.5 191.4 191.7 191.9 191.9 192.2 192.7 193.0 193.1 193.2
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 184.0 183.9 183.9 183.9 183.9 183.8 183.8 183.8 183.8 183.7 183.7
## 52% 186.5 186.5 186.4 186.5 186.3 186.4 186.4 186.4 186.4 186.5 186.4
## 60% 187.5 187.4 187.3 187.3 187.3 187.3 187.4 187.4 187.5 187.4 187.4
## 70% 188.4 188.4 188.6 188.6 188.7 188.8 188.6 188.7 188.7 188.9 188.9
## 95% 193.4 193.3 193.5 194.0 193.9 194.0 194.1 194.2 194.7 194.7 194.6
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 183.6 183.5 183.4 183.3 183.4 183.3 183.2 183.1 183.0 183.0 183.1
## 52% 186.5 186.5 186.4 186.5 186.4 186.4 186.2 186.5 186.5 186.5 186.3
## 60% 187.5 187.7 187.5 187.4 187.5 187.6 187.5 187.7 187.7 187.8 187.7
## 70% 188.8 188.9 189.0 189.2 189.1 189.1 189.1 189.3 189.4 189.5 189.3
## 95% 194.6 195.1 195.3 195.3 195.6 195.7 195.7 195.9 196.2 196.0 196.3
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 183.0 183.0 183.0 183.0 182.9 183.0 182.9 182.7 182.7 182.8 182.9
## 52% 186.3 186.3 186.4 186.3 186.4 186.1 186.4 186.4 186.5 186.6 186.4
## 60% 187.7 187.7 187.7 187.7 187.8 187.6 188.0 187.8 187.9 187.8 187.9
## 70% 189.5 189.5 189.5 189.6 189.5 189.4 189.6 189.8 189.8 189.8 189.8
## 95% 195.6 196.0 195.9 196.0 196.4 196.3 196.3 196.6 196.8 197.2 197.1
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 182.9 182.9 183.0 183.0 182.9 182.8 182.7 182.9 182.8 182.9 182.7
## 52% 186.4 186.5 186.1 186.5 186.6 186.4 186.6 186.5 186.4 186.2 186.2
## 60% 187.9 187.9 187.8 187.8 187.8 187.8 187.9 187.9 187.7 187.9 187.7
## 70% 189.7 189.7 189.9 189.9 189.9 189.8 189.8 189.9 189.8 189.9 190.0
## 95% 197.1 197.0 197.0 197.1 197.3 197.2 197.6 197.7 197.9 198.3 198.3
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 182.7 182.6 182.5 182.7 182.5 182.4 182.3 182.2 182.3 182.2 182.3
## 52% 186.3 186.3 186.4 186.1 186.4 186.6 186.5 186.4 186.7 186.6 186.8
## 60% 187.8 187.8 187.8 187.7 187.8 188.0 188.0 188.1 188.3 188.3 188.5
## 70% 190.0 190.0 190.1 190.2 190.1 190.3 190.1 190.1 190.2 190.4 190.1
## 95% 198.2 198.3 198.7 198.6 198.8 198.9 199.2 198.9 199.1 199.2 199.4
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 182.3 182.4 182.2 182.3 182.1 182.1 182.1 182.0 181.9 181.6 181.8
## 52% 186.7 186.6 186.6 186.6 186.5 186.6 186.6 186.5 186.7 186.6 186.7
## 60% 188.2 188.4 188.3 188.2 188.2 188.1 188.1 188.2 188.4 188.2 188.1
## 70% 190.2 190.3 190.5 190.4 190.4 190.3 190.6 190.7 190.5 190.5 190.5
## 95% 199.3 199.6 199.5 199.9 199.8 200.2 200.6 200.5 200.6 200.6 200.8
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 181.8 181.6 181.7 181.6 181.5 181.6 181.6 181.6 181.5 181.4 181.5
## 52% 186.7 187.1 186.9 186.9 186.9 186.9 186.9 186.8 186.7 186.4 186.5
## 60% 188.2 188.3 188.4 188.2 188.2 188.3 188.3 188.4 188.3 188.3 188.3
## 70% 190.5 190.5 190.4 190.4 190.6 190.6 190.4 190.2 190.3 190.4 190.7
## 95% 201.1 200.8 200.7 200.7 200.5 200.8 200.6 200.7 201.1 201.2 201.0
## [,100]
## 31% 181.2
## 52% 186.6
## 60% 188.3
## 70% 190.4
## 95% 201.1
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.8 187.7 187.6 187.4 187.4 187.3 187.1 186.9
## 52% 188.8 188.9 188.9 188.9 188.9 189.0 189.1 189.1 189.1 189.1 189.0
## 60% 188.8 189.1 189.2 189.3 189.4 189.4 189.6 189.7 189.7 189.8 189.8
## 70% 188.8 189.4 189.7 189.8 190.0 190.2 190.3 190.5 190.6 190.8 190.7
## 95% 188.8 190.8 191.5 192.1 192.5 192.8 193.4 193.6 193.9 194.3 194.7
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 187.0 186.9 186.8 186.7 186.6 186.4 186.4 186.3 186.2 186.2 186.0
## 52% 188.9 189.1 189.0 189.0 188.9 189.0 189.1 189.0 189.0 189.0 188.9
## 60% 189.7 189.8 189.8 189.8 190.0 190.1 189.9 189.9 189.9 190.0 189.9
## 70% 190.9 190.9 190.9 190.9 191.2 191.2 191.3 191.4 191.5 191.4 191.4
## 95% 194.7 195.5 195.7 196.1 196.3 196.4 196.7 197.1 197.1 197.3 197.8
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 186.1 185.9 186.0 185.9 186.1 185.9 185.9 185.9 185.7 185.7 185.6
## 52% 188.9 189.0 188.9 189.2 189.1 189.3 189.3 189.2 189.1 188.9 189.0
## 60% 190.0 190.0 190.1 190.2 190.3 190.3 190.3 190.4 190.2 190.2 190.4
## 70% 191.5 191.5 191.6 191.6 191.7 191.8 191.7 191.8 192.0 192.2 192.2
## 95% 197.9 197.9 198.3 198.2 198.4 198.4 198.3 198.5 198.8 198.9 198.8
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 185.4 185.3 185.2 185.2 185.3 185.2 185.3 185.1 185.0 185.1 185.0
## 52% 189.1 189.1 189.2 189.2 189.1 189.0 189.0 189.0 189.0 189.1 189.1
## 60% 190.4 190.4 190.3 190.4 190.4 190.3 190.4 190.7 190.5 190.5 190.4
## 70% 192.3 192.3 192.3 192.3 192.5 192.2 192.4 192.6 192.8 192.6 192.7
## 95% 198.9 198.9 199.1 199.6 199.8 200.2 200.5 200.7 200.6 200.9 201.2
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 185.1 185.0 184.9 184.9 184.8 184.8 184.9 184.9 184.7 184.8 184.5
## 52% 188.9 189.0 189.2 189.2 189.3 189.2 189.3 189.4 189.3 189.1 189.2
## 60% 190.6 190.5 190.6 190.9 190.9 190.9 190.8 190.9 190.9 190.8 191.0
## 70% 192.6 192.6 192.8 193.0 192.8 193.0 192.9 192.9 192.9 193.0 193.1
## 95% 200.9 201.3 201.7 202.1 201.7 201.9 202.2 202.3 202.4 202.7 202.8
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 184.4 184.5 184.3 184.5 184.4 184.4 184.5 184.5 184.3 184.6 184.4
## 52% 189.3 189.1 189.1 188.9 189.1 189.2 188.8 189.0 189.2 189.1 189.2
## 60% 190.9 190.7 190.7 190.7 190.8 190.7 190.5 190.8 190.8 191.0 190.7
## 70% 193.2 193.3 193.4 193.3 193.2 193.4 193.3 193.3 193.3 193.4 193.5
## 95% 202.9 203.0 203.1 203.2 203.1 203.0 203.1 203.2 203.2 203.6 203.6
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 184.3 184.3 184.1 184.2 184.0 183.8 183.8 183.8 183.8 184.1 184.1
## 52% 189.2 189.2 189.2 189.4 189.3 189.4 189.3 189.2 189.2 189.2 189.2
## 60% 190.9 190.8 190.9 191.2 191.3 191.1 191.2 191.3 191.3 191.1 191.3
## 70% 193.7 193.4 193.4 193.3 193.9 193.7 193.6 193.6 193.7 193.6 194.0
## 95% 203.6 203.9 204.3 204.4 204.4 203.9 204.1 204.1 204.5 204.4 204.6
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 184.0 183.8 183.9 183.7 183.8 183.5 183.5 183.7 183.4 183.3 183.2
## 52% 189.2 189.2 189.2 189.2 189.1 189.3 189.2 188.9 189.0 189.1 188.8
## 60% 191.3 191.1 191.2 191.0 191.2 191.1 191.3 191.1 191.4 191.5 191.2
## 70% 193.8 193.7 193.8 193.8 193.9 193.5 193.8 193.8 193.9 194.1 193.9
## 95% 205.0 204.4 204.8 204.8 204.5 204.7 205.1 205.0 205.3 205.2 205.4
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 183.2 183.2 183.1 183.2 183.1 183.0 183.1 183.3 183.5 183.2 183.2
## 52% 188.8 188.8 188.8 188.8 188.6 188.6 188.8 188.8 188.9 188.8 188.8
## 60% 191.0 191.1 191.0 190.9 190.8 190.5 190.8 190.8 190.7 190.8 191.1
## 70% 193.8 193.7 193.7 193.4 193.4 193.6 193.4 193.6 193.6 193.5 193.6
## 95% 205.5 205.6 205.4 205.5 206.0 206.3 206.0 206.3 206.0 206.1 206.1
## [,100]
## 31% 183.3
## 52% 188.9
## 60% 191.2
## 70% 193.8
## 95% 206.2
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.8 506.6 506.3 506.0 505.9 505.8 505.4 505.3 505.2 504.9
## 52% 507.7 507.8 507.9 507.9 508.0 508.1 508.1 508.3 508.2 508.1 508.3
## 60% 507.7 508.2 508.4 508.6 508.7 508.9 509.0 509.1 509.3 509.2 509.4
## 70% 507.7 508.7 509.1 509.4 509.7 509.9 510.1 510.4 510.6 510.8 510.7
## 95% 507.7 510.7 511.7 513.0 513.9 514.5 515.1 515.8 516.2 517.0 517.2
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 505.0 504.8 504.7 504.6 504.6 504.4 504.4 504.4 503.9 504.0 503.6
## 52% 508.2 508.2 508.1 508.0 508.0 508.0 508.0 508.2 508.2 508.1 508.1
## 60% 509.4 509.3 509.4 509.6 509.5 509.8 509.9 509.9 509.7 509.9 509.9
## 70% 511.0 510.9 511.3 511.5 511.5 511.8 511.8 511.9 511.9 512.0 512.0
## 95% 517.5 518.2 518.7 518.8 519.3 520.2 520.3 520.8 521.4 521.5 522.0
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 503.5 503.5 503.4 503.4 503.3 503.1 503.2 503.1 503.1 502.5 502.8
## 52% 508.2 508.4 508.2 508.2 508.3 508.2 508.3 508.4 508.2 508.2 508.0
## 60% 510.0 509.9 510.0 510.0 510.2 510.2 509.9 510.2 510.3 510.2 510.1
## 70% 512.2 512.4 512.4 512.3 512.7 512.5 512.4 512.7 512.8 512.6 513.1
## 95% 522.4 523.0 522.8 522.7 522.5 522.8 523.6 523.8 524.1 524.9 525.3
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 502.6 502.6 502.7 502.6 502.4 501.9 501.9 502.3 502.0 502.1 502.0
## 52% 508.1 508.2 508.0 507.7 508.1 508.0 507.9 508.0 508.0 507.8 508.0
## 60% 510.1 510.0 510.0 510.1 510.2 510.2 510.4 510.6 510.4 510.5 510.5
## 70% 512.9 513.0 513.2 513.1 513.3 513.5 513.6 513.9 513.6 513.6 513.6
## 95% 525.4 525.7 526.9 525.9 525.6 525.9 526.2 526.4 526.6 527.6 528.1
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 502.0 501.4 501.3 501.2 501.1 500.9 501.3 501.2 501.1 501.1 501.0
## 52% 508.1 508.4 508.4 508.5 508.1 507.8 507.8 507.8 507.8 507.6 507.9
## 60% 510.6 510.5 510.8 510.7 510.7 510.2 510.3 510.3 510.3 510.1 510.2
## 70% 513.3 513.8 513.7 514.1 514.2 514.1 514.2 514.1 514.3 513.8 513.9
## 95% 527.1 527.5 527.3 527.9 529.0 529.2 529.3 530.2 529.8 530.0 529.8
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 500.8 500.9 500.5 500.6 500.5 500.5 500.5 500.2 500.3 500.2 500.0
## 52% 507.7 507.5 507.7 507.8 507.6 507.4 507.9 507.7 507.9 507.7 508.1
## 60% 510.4 510.1 510.5 510.6 510.6 510.6 510.3 510.5 510.8 510.3 510.6
## 70% 514.1 514.0 514.4 514.5 514.3 514.2 514.3 514.2 514.5 514.3 514.6
## 95% 530.6 530.9 531.1 532.0 532.2 532.1 532.9 532.8 533.2 534.0 534.3
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 499.7 499.5 499.4 499.5 499.2 498.9 499.4 499.4 499.2 499.1 498.6
## 52% 507.8 508.0 508.0 508.4 508.6 508.5 508.4 508.0 508.1 507.8 507.9
## 60% 510.9 511.2 511.2 510.8 511.1 510.9 511.2 511.4 511.1 511.2 511.5
## 70% 514.6 514.8 515.2 515.2 515.3 515.5 515.3 515.5 515.3 516.2 516.1
## 95% 534.1 534.6 534.7 534.1 534.2 534.9 534.9 534.8 535.4 535.6 535.8
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 499.2 498.7 498.5 498.6 498.5 498.3 498.4 498.5 498.2 498.3 498.5
## 52% 507.9 508.3 508.1 508.3 508.4 507.8 507.7 508.1 508.4 508.0 507.9
## 60% 511.2 511.4 510.9 511.4 511.4 512.0 512.0 511.9 512.0 511.4 511.1
## 70% 515.5 515.3 515.4 515.4 515.8 515.6 515.7 515.9 515.8 516.0 515.7
## 95% 534.9 534.9 535.9 536.5 536.2 536.8 536.2 536.1 536.4 536.7 536.3
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 498.2 498.4 498.0 497.9 497.8 498.1 498.1 498.5 498.3 498.7 498.4
## 52% 507.6 507.7 508.0 508.4 508.2 508.3 508.2 508.0 508.0 508.1 508.1
## 60% 511.1 511.7 511.8 511.5 511.7 511.6 511.3 511.8 511.8 512.1 511.6
## 70% 516.2 515.6 515.7 516.0 516.2 516.2 516.2 516.7 516.6 516.7 517.0
## 95% 536.5 536.1 536.4 537.3 537.9 537.8 537.1 537.3 537.4 538.4 537.8
## [,100]
## 31% 498.7
## 52% 508.3
## 60% 511.3
## 70% 517.2
## 95% 537.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(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.7 278.4 278.2 278.0 277.8 277.6 277.3 277.3 277.1 277.0
## 52% 279.5 279.6 279.6 279.6 279.4 279.5 279.6 279.7 279.5 279.5 279.6
## 60% 279.5 279.9 280.0 280.0 280.0 280.1 280.2 280.4 280.4 280.4 280.4
## 70% 279.5 280.2 280.5 280.6 280.8 281.0 281.1 281.5 281.4 281.5 281.6
## 95% 279.5 281.8 282.8 283.2 283.8 284.5 284.9 285.3 285.6 286.0 286.5
## [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
## 31% 276.8 276.7 276.5 276.5 276.5 276.2 276.2 276.0 276.0 275.9 275.8
## 52% 279.6 279.5 279.4 279.4 279.4 279.5 279.5 279.4 279.4 279.2 279.2
## 60% 280.5 280.4 280.5 280.4 280.3 280.6 280.5 280.7 280.6 280.5 280.5
## 70% 281.7 281.8 281.8 281.7 281.7 282.0 282.2 282.1 282.2 282.2 282.3
## 95% 286.5 287.2 287.4 287.9 288.2 288.5 288.8 288.9 289.0 289.8 289.8
## [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33]
## 31% 275.8 275.7 275.7 275.8 275.5 275.8 275.7 275.4 275.1 274.8 274.7
## 52% 279.3 279.4 279.5 279.4 279.5 279.3 279.2 279.2 279.0 279.2 279.2
## 60% 280.7 280.8 280.9 281.0 281.0 281.0 281.0 280.9 280.8 280.7 280.7
## 70% 282.4 282.5 282.6 282.7 283.0 283.3 283.3 283.2 283.1 282.9 283.1
## 95% 290.0 290.1 291.2 291.0 291.0 291.0 291.3 291.5 292.2 292.6 292.9
## [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44]
## 31% 274.7 274.7 274.7 274.6 274.5 274.3 274.2 274.3 274.4 274.0 274.4
## 52% 279.2 279.2 279.4 279.4 279.4 279.3 279.2 279.4 279.3 279.5 279.5
## 60% 281.0 281.1 281.0 281.1 281.2 281.0 281.1 281.4 281.4 281.2 281.6
## 70% 283.5 283.5 283.5 283.6 283.7 283.9 284.0 283.9 284.0 283.9 283.9
## 95% 292.8 293.3 293.2 293.2 293.5 293.2 294.1 294.2 294.4 294.8 294.7
## [,45] [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55]
## 31% 274.4 274.1 274.2 274.5 274.1 274.2 274.0 274.0 273.9 273.8 273.9
## 52% 279.4 279.5 279.6 279.5 279.4 279.5 279.5 279.5 279.4 279.7 279.2
## 60% 281.3 281.5 281.5 281.6 281.3 281.5 281.6 281.6 281.6 281.6 281.6
## 70% 284.2 284.2 284.2 284.3 284.2 284.3 284.4 284.7 284.4 284.0 284.2
## 95% 294.7 294.4 294.7 294.5 294.6 295.0 295.2 295.5 295.7 296.3 296.5
## [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66]
## 31% 273.9 273.6 273.8 273.5 273.6 273.1 273.2 273.2 273.5 273.3 273.2
## 52% 279.3 279.2 279.1 279.0 278.9 278.7 279.1 279.3 279.1 279.1 279.4
## 60% 281.5 281.3 281.3 281.4 281.4 281.6 281.4 281.3 281.3 281.4 281.2
## 70% 284.2 284.2 284.1 284.3 284.6 284.3 284.6 284.8 284.6 284.6 284.7
## 95% 296.8 296.8 297.1 297.3 297.9 297.6 298.1 297.8 298.5 298.3 298.7
## [,67] [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77]
## 31% 273.0 273.2 273.0 273.5 273.2 273.2 272.9 273.0 273.1 272.7 272.8
## 52% 279.1 279.2 279.4 279.3 279.2 279.3 279.3 279.4 279.5 279.4 279.1
## 60% 281.6 281.4 281.7 281.8 281.8 281.8 281.6 281.7 281.7 281.7 281.8
## 70% 284.7 285.0 285.1 285.0 285.2 285.3 285.0 285.1 285.0 285.2 285.3
## 95% 298.4 298.7 298.5 298.9 298.8 298.8 299.5 299.3 299.5 299.8 299.7
## [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88]
## 31% 272.6 272.5 272.7 272.4 272.5 272.2 272.3 272.2 271.7 271.9 271.9
## 52% 279.3 279.5 279.6 279.7 279.6 279.5 279.4 279.9 279.8 279.7 279.5
## 60% 281.7 281.9 281.9 281.8 282.2 282.3 282.6 282.6 282.3 282.3 282.4
## 70% 285.6 285.9 285.6 285.5 285.6 285.7 286.0 286.1 286.2 286.3 286.3
## 95% 299.3 299.9 299.8 300.0 300.0 300.8 300.4 300.4 300.9 301.0 300.7
## [,89] [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99]
## 31% 271.8 272.3 272.3 272.5 272.4 272.1 272.0 272.1 272.4 272.3 272.1
## 52% 279.9 279.6 279.6 279.9 280.1 280.1 280.1 279.9 279.8 280.1 280.0
## 60% 282.4 282.3 282.5 282.6 282.7 282.7 283.0 283.1 282.8 283.1 282.9
## 70% 286.0 285.8 286.1 285.9 286.1 286.2 286.2 286.0 286.2 286.6 286.6
## 95% 301.1 301.2 301.3 301.2 301.7 301.3 302.0 301.5 302.0 302.1 303.0
## [,100]
## 31% 272.3
## 52% 280.0
## 60% 282.9
## 70% 286.3
## 95% 302.4
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.