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); Second, we manually use the function from “Calibrating the Ornstein-Uhlenbeck (Vasicek) model”http://www.sitmo.com/article/calibrating-the-ornstein-uhlenbeck-model/
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)
## Loading required package: stats4
require(sde)
## Loading required package: sde
## Loading required package: MASS
## Loading required package: fda
## Loading required package: splines
## Loading required package: Matrix
##
## Attaching package: 'fda'
##
## The following object is masked from 'package:graphics':
##
## matplot
##
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## To check the errata corrige of the book, type vignette("sde.errata")
# 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)
## [1] 169.9 170.7 174.5 167.2 165.2 167.4
X <- ts(data = wheat)
fitwheat <- mle(OU.lik, start = list(theta1 = 20, theta2 = 0.1, theta3 = 10),
method = "L-BFGS-B", lower = c(-Inf, 0, 0))
## Error: non-numeric argument to binary operator
summary(fitwheat)
## Error: error in evaluating the argument 'object' in selecting a method for function 'summary': Error: object 'fitwheat' not found
coef(fitwheat)[[1]]/coef(fitwheat)[[2]]
## Error: error in evaluating the argument 'object' in selecting a method for function 'coef': Error: object 'fitwheat' not found
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
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 \]
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: expm
##
## 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 = 220)
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] [,12]
## 31% 220 218.9 218.4 218.2 217.9 217.8 217.5 217.7 217.3 216.8 216.6 216.7
## 52% 220 220.2 220.1 220.1 220.2 220.3 220.4 220.2 220.5 220.4 220.5 220.8
## 60% 220 220.6 220.8 221.0 221.2 221.3 221.5 221.3 221.6 221.6 222.0 222.2
## 70% 220 221.3 221.7 221.9 222.3 222.8 222.9 223.1 223.2 223.2 223.7 224.1
## 95% 220 223.4 225.0 226.1 227.2 227.6 229.0 229.5 230.4 231.5 231.9 231.8
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## 31% 216.4 216.3 216.4 216.1 215.9 216.1 215.8 215.9 216.0 216.0 215.8
## 52% 220.7 220.7 220.6 220.9 221.1 221.2 221.0 221.1 220.9 221.0 221.2
## 60% 222.2 222.0 222.3 222.7 223.0 223.1 223.1 222.9 222.9 223.2 223.7
## 70% 223.9 224.2 224.8 225.0 225.1 225.6 225.6 225.6 226.0 226.0 226.2
## 95% 232.7 233.0 233.6 233.6 235.0 234.9 235.5 235.6 235.9 236.7 236.9
## [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## 31% 216.2 215.9 215.9 215.7 215.9 216.2 215.8 215.9 215.4 215.8 215.3
## 52% 221.6 221.6 221.5 221.7 221.8 221.9 221.9 222.1 222.2 222.4 222.7
## 60% 223.5 223.9 224.0 223.8 224.0 224.2 224.2 224.7 224.9 224.8 224.9
## 70% 226.1 226.7 226.8 227.3 227.0 227.1 227.2 227.7 227.7 228.1 228.1
## 95% 237.2 238.3 238.4 238.7 239.0 239.7 240.4 240.2 240.6 241.2 241.9
## [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
## 31% 215.5 215.7 215.3 215.3 215.2 214.9 215.1 214.9 214.8 214.5 214.4
## 52% 222.4 222.7 222.7 222.5 222.8 223.0 222.8 222.8 222.5 222.7 222.9
## 60% 225.3 225.1 225.4 225.5 225.7 225.6 225.8 226.0 225.9 225.8 225.9
## 70% 228.0 228.0 228.2 228.6 229.1 229.4 229.7 229.5 229.7 229.5 229.7
## 95% 241.7 241.8 241.6 242.1 242.4 243.1 243.9 244.2 244.1 245.2 245.2
## [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## 31% 214.8 214.8 214.5 214.3 214.3 214.7 214.8 214.9 214.5 214.6 214.5
## 52% 222.6 222.7 222.9 222.4 222.7 222.6 222.7 223.1 223.6 223.2 223.2
## 60% 225.8 225.8 225.8 225.9 225.7 226.0 226.2 226.1 226.3 226.5 226.8
## 70% 229.6 229.6 229.3 229.6 229.2 229.7 229.9 230.1 230.0 230.5 230.6
## 95% 245.4 244.8 245.6 246.3 246.2 245.9 245.9 246.7 247.1 246.9 247.7
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67]
## 31% 214.7 214.2 214.3 214.8 214.7 214.4 214.8 214.7 214.0 213.9 214.1
## 52% 222.9 223.2 223.5 223.6 223.3 223.2 223.5 223.4 223.8 223.8 223.9
## 60% 226.4 226.3 226.2 226.9 227.0 227.1 226.9 227.0 226.7 226.9 227.1
## 70% 230.6 230.3 230.7 230.9 231.0 231.2 231.1 231.5 232.1 232.0 231.9
## 95% 247.9 247.6 247.5 248.5 249.1 249.0 249.1 249.2 249.7 250.1 250.5
## [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78]
## 31% 214.1 214.3 213.9 213.4 213.8 213.7 214.1 214.1 214.1 213.9 213.6
## 52% 223.5 224.0 224.0 223.8 224.0 224.0 223.8 224.2 224.4 224.6 224.5
## 60% 226.9 227.0 227.2 227.4 227.6 227.3 227.7 227.7 227.7 227.8 228.0
## 70% 232.0 232.3 232.2 232.7 232.4 232.5 232.5 232.9 232.8 232.8 233.1
## 95% 250.6 250.5 251.3 251.5 251.5 251.9 251.6 251.8 252.0 252.7 252.8
## [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89]
## 31% 213.7 214.1 214.1 213.8 213.6 213.5 213.5 213.2 213.6 213.5 213.4
## 52% 224.6 224.7 224.7 224.6 224.4 224.5 224.5 224.3 223.9 223.8 223.9
## 60% 228.2 228.3 228.3 228.3 228.2 228.3 228.2 228.3 228.6 227.7 228.1
## 70% 233.2 233.0 232.9 233.3 233.3 233.1 233.1 233.5 233.8 233.7 234.1
## 95% 252.9 253.2 254.0 253.2 254.2 254.0 254.5 254.6 255.7 256.4 256.2
## [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## 31% 213.1 213.4 213.1 212.7 213.3 213.0 213.1 212.9 213.0 212.8 213.1
## 52% 224.5 224.5 224.4 224.8 224.9 224.6 225.0 225.2 225.3 225.7 225.0
## 60% 228.2 228.3 229.0 228.7 228.8 228.5 228.7 228.6 228.8 228.8 229.3
## 70% 233.8 233.4 233.9 234.0 234.3 234.2 234.1 234.2 234.5 234.5 234.7
## 95% 256.1 257.1 256.3 257.0 257.0 257.8 258.5 258.8 259.1 259.9 259.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 = 188)
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] [,12]
## 31% 188 187.6 187.3 187.0 186.9 186.8 186.8 186.6 186.6 186.5 186.3 186.2
## 52% 188 188.0 188.0 188.0 188.0 188.1 188.0 188.0 188.0 187.9 188.1 188.1
## 60% 188 188.2 188.3 188.3 188.4 188.5 188.5 188.5 188.5 188.7 188.7 188.7
## 70% 188 188.5 188.6 188.8 189.0 189.0 189.1 189.2 189.2 189.4 189.5 189.5
## 95% 188 189.6 190.0 190.6 191.1 191.2 191.5 192.1 192.1 192.4 192.8 192.8
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## 31% 186.2 186.2 186.0 185.9 185.8 185.9 185.8 185.8 185.8 185.6 185.8
## 52% 188.1 188.0 188.1 188.0 188.0 188.0 188.0 188.0 188.1 188.2 188.1
## 60% 188.8 188.8 188.8 188.9 188.8 188.8 188.9 188.9 189.1 189.0 189.0
## 70% 189.5 189.6 189.8 189.7 189.7 189.9 189.9 189.9 190.1 190.2 190.1
## 95% 193.3 193.4 193.6 193.7 193.6 194.2 194.5 194.6 194.9 194.6 195.0
## [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## 31% 185.7 185.6 185.5 185.6 185.5 185.3 185.3 185.3 185.1 185.0 185.0
## 52% 188.1 188.0 187.9 187.9 187.9 188.0 188.1 188.0 188.2 188.1 188.0
## 60% 188.9 188.9 188.8 189.1 188.9 189.0 189.0 189.2 189.1 189.0 189.2
## 70% 190.2 190.2 190.3 190.4 190.3 190.2 190.5 190.4 190.4 190.3 190.5
## 95% 195.4 195.4 195.3 195.5 195.5 195.3 195.7 195.7 196.0 196.2 196.3
## [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
## 31% 185.0 185.0 184.8 184.8 184.7 184.6 184.8 184.7 184.5 184.4 184.5
## 52% 187.9 187.8 187.9 187.9 187.7 187.9 188.0 187.9 187.9 187.8 187.9
## 60% 189.0 189.0 189.1 189.1 189.2 189.1 189.1 189.3 189.1 189.3 189.1
## 70% 190.6 190.7 190.4 190.8 190.8 190.4 190.7 190.6 190.5 190.5 190.7
## 95% 196.7 196.7 196.9 197.2 197.0 197.3 197.4 197.6 197.8 197.9 197.8
## [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## 31% 184.6 184.6 184.5 184.4 184.5 184.4 184.3 184.4 184.6 184.4 184.4
## 52% 188.0 187.9 188.0 188.0 188.0 187.9 187.9 187.9 187.8 187.9 188.0
## 60% 189.2 189.1 189.1 189.3 189.2 189.2 189.1 189.1 189.2 189.3 189.2
## 70% 190.7 190.7 190.8 190.8 190.9 190.8 190.9 190.8 191.0 191.0 191.0
## 95% 198.2 198.0 198.0 198.6 198.8 198.4 198.8 198.6 198.9 199.0 198.9
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67]
## 31% 184.4 184.3 184.3 184.3 184.3 184.3 184.2 184.2 184.2 184.2 184.1
## 52% 187.9 187.8 187.8 187.8 187.7 187.8 187.8 187.8 187.9 187.9 188.0
## 60% 189.2 189.1 189.2 189.1 189.2 189.2 189.3 189.4 189.3 189.3 189.2
## 70% 191.2 191.0 191.0 191.1 191.2 191.2 191.2 191.3 191.2 191.2 191.2
## 95% 198.9 199.4 199.3 199.2 199.4 199.7 199.7 199.5 199.8 199.5 199.6
## [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78]
## 31% 183.9 183.9 183.8 183.8 183.7 183.7 183.5 183.8 183.5 183.7 183.8
## 52% 188.0 187.8 187.8 187.7 188.0 187.8 187.9 187.9 187.9 187.8 188.0
## 60% 189.3 189.3 189.3 189.5 189.4 189.5 189.5 189.5 189.5 189.3 189.4
## 70% 191.1 191.4 191.3 191.4 191.4 191.4 191.4 191.3 191.6 191.5 191.7
## 95% 199.7 199.9 199.8 199.8 200.1 200.2 200.7 200.7 200.9 201.0 200.9
## [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89]
## 31% 183.5 183.5 183.4 183.4 183.4 183.6 183.4 183.6 183.6 183.7 183.6
## 52% 188.0 188.2 188.2 188.0 188.2 188.2 188.2 188.2 188.3 188.4 188.7
## 60% 189.5 189.6 189.5 189.8 189.8 190.0 189.7 189.8 189.8 189.8 190.0
## 70% 191.8 191.5 191.7 191.7 191.6 192.1 192.0 192.0 192.0 192.0 192.3
## 95% 201.2 200.9 200.9 201.2 201.3 201.2 201.3 201.5 201.6 201.7 201.9
## [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## 31% 183.6 183.6 183.5 183.6 183.4 183.6 183.3 183.5 183.4 183.5 183.1
## 52% 188.5 188.5 188.7 188.4 188.5 188.2 188.4 188.3 188.5 188.4 188.4
## 60% 190.0 190.1 190.0 190.0 190.0 190.0 190.1 190.3 190.0 190.0 190.0
## 70% 192.3 192.3 192.4 192.5 192.4 192.4 192.5 192.5 192.4 192.4 192.7
## 95% 202.3 202.1 202.1 202.2 202.3 203.1 203.2 202.9 203.0 203.4 203.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 = 188)
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] [,12]
## 31% 188 187.4 187.2 186.9 187.0 186.7 186.6 186.4 186.4 186.3 186.2 186.2
## 52% 188 188.1 188.1 188.0 188.1 188.1 188.0 188.0 188.0 188.0 188.2 188.2
## 60% 188 188.3 188.5 188.5 188.6 188.6 188.6 188.6 188.6 188.6 188.9 188.9
## 70% 188 188.6 188.9 189.0 189.2 189.3 189.3 189.3 189.4 189.6 189.7 189.9
## 95% 188 189.9 190.6 191.1 191.8 192.1 192.3 192.8 193.1 193.4 193.7 194.1
## [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23]
## 31% 186.2 186.2 186.0 186.0 186.0 185.9 185.9 185.8 185.8 185.9 185.6
## 52% 188.2 188.3 188.5 188.3 188.3 188.4 188.4 188.4 188.4 188.4 188.4
## 60% 189.0 189.0 189.1 189.2 189.3 189.4 189.5 189.5 189.4 189.4 189.4
## 70% 190.0 190.2 190.4 190.5 190.5 190.6 190.8 190.8 191.0 191.2 191.1
## 95% 194.6 194.8 195.1 195.4 195.9 195.9 196.2 196.5 196.8 196.9 197.0
## [,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## 31% 185.7 185.5 185.3 185.4 185.3 185.4 185.2 185.2 185.2 185.0 185.1
## 52% 188.4 188.5 188.4 188.5 188.5 188.3 188.3 188.4 188.3 188.5 188.4
## 60% 189.5 189.6 189.6 189.6 189.7 189.5 189.4 189.6 189.7 189.6 189.8
## 70% 191.0 191.0 191.1 191.3 191.4 191.3 191.2 191.6 191.4 191.4 191.5
## 95% 197.0 197.1 197.7 197.9 197.8 197.9 198.1 198.5 198.4 198.5 198.4
## [,35] [,36] [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
## 31% 185.0 185.0 185.0 184.9 184.8 184.7 184.5 184.4 184.5 184.6 184.5
## 52% 188.4 188.5 188.5 188.5 188.5 188.6 188.5 188.7 188.5 188.4 188.4
## 60% 189.6 189.8 189.7 189.8 189.8 190.0 189.9 190.0 189.8 190.0 189.9
## 70% 191.3 191.5 191.6 191.6 191.6 191.5 191.6 191.8 191.9 191.8 191.7
## 95% 198.7 199.0 199.3 199.4 199.4 199.2 199.2 199.7 199.3 199.5 199.8
## [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54] [,55] [,56]
## 31% 184.8 184.6 184.4 184.5 184.4 184.3 184.2 184.1 184.3 184.3 184.1
## 52% 188.5 188.4 188.4 188.5 188.4 188.5 188.5 188.2 188.4 188.6 188.9
## 60% 189.8 189.9 189.8 190.1 190.0 189.8 189.8 189.8 190.0 190.2 190.3
## 70% 191.9 192.0 192.2 192.1 192.2 192.1 192.0 192.2 192.1 192.2 192.4
## 95% 200.3 200.3 200.5 200.7 200.9 200.6 200.7 201.0 201.4 201.1 201.1
## [,57] [,58] [,59] [,60] [,61] [,62] [,63] [,64] [,65] [,66] [,67]
## 31% 184.1 183.9 183.8 183.9 183.8 183.7 183.8 183.6 183.7 183.9 183.9
## 52% 188.6 188.8 188.9 188.7 188.7 188.8 188.8 188.7 188.8 188.8 188.7
## 60% 190.5 190.8 190.8 190.8 190.9 190.7 190.5 190.7 190.4 190.5 190.7
## 70% 192.6 192.7 192.9 192.8 193.0 192.9 192.8 193.0 193.3 193.4 193.3
## 95% 201.2 202.1 201.6 202.3 202.4 202.8 202.8 203.0 202.7 203.3 203.7
## [,68] [,69] [,70] [,71] [,72] [,73] [,74] [,75] [,76] [,77] [,78]
## 31% 183.8 183.7 183.5 183.7 183.5 183.5 183.6 183.5 183.4 183.3 183.3
## 52% 188.9 188.9 189.1 189.1 189.1 189.0 189.1 189.3 189.2 189.1 189.1
## 60% 190.8 190.6 190.8 190.9 191.0 191.0 191.0 191.0 191.3 191.2 191.0
## 70% 193.2 193.5 193.5 193.5 193.4 193.5 193.6 193.7 193.9 193.9 194.0
## 95% 203.6 203.7 203.8 204.0 203.8 204.2 204.0 203.9 203.8 204.6 203.9
## [,79] [,80] [,81] [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89]
## 31% 183.3 183.1 183.2 183.2 183.0 182.9 183.0 182.8 182.7 182.8 182.7
## 52% 189.0 189.0 188.9 189.4 189.3 189.2 189.2 189.2 188.9 189.3 189.1
## 60% 191.2 191.1 191.3 191.4 191.4 191.5 191.3 191.3 191.5 191.5 191.4
## 70% 194.0 194.3 194.0 194.2 194.1 194.1 194.4 194.2 194.4 194.3 194.5
## 95% 204.4 204.6 205.1 204.8 205.5 205.3 205.2 205.4 205.4 205.8 206.2
## [,90] [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## 31% 182.6 182.8 182.7 182.8 182.5 182.8 182.6 182.7 182.6 182.7 182.6
## 52% 189.2 189.2 189.0 189.3 189.4 189.4 189.5 189.3 189.3 189.3 189.4
## 60% 191.2 191.3 191.3 191.1 191.4 191.5 191.4 191.5 191.6 191.8 191.5
## 70% 194.3 194.3 194.3 194.4 194.3 194.3 194.3 194.2 194.6 194.4 194.5
## 95% 206.1 205.8 206.1 206.3 206.5 206.4 206.8 206.7 207.1 207.6 207.7
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)
##
## T set to = 1000.000000
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)
##
## T set to = 1000.000000
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)
##
## T set to = 1000.000000
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)
##
## T set to = 1000.000000
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)
##
## T set to = 1000.000000
plot(Xdry.peas, main = " Ornstein - Uhlenbeck ")
