obtenemos los datos desde el sitio de yahoo

SBUXdata <- pdfetch_YAHOO("SBUX",from = c("2019-01-01"),to = c("2020-01-01"), interval = '1d')

Obtenemos la columna que es de nuestro interés

Novocure <- SBUXdata[,4]
length(Novocure)
## [1] 252

Convertimos en una serie temporal

tsNovoCure <- ts(Novocure, start = c(2019,1),frequency=365)

Se calculan las diferencias de la serie de datos con logaritmo

#Calculo parámetros iniciales manera 1

l_NovoCure<-diff(log(tsNovoCure))
plot(l_NovoCure)#Calculo parámetros iniciales manera 1

Delta <- 1/365
alpha <- mean(l_NovoCure)/Delta
sigma <- sqrt(var(l_NovoCure)/Delta)
mu <- alpha +0.5*sigma^2
x0<-tsNovoCure[1]
x0
## [1] 64.32
sigma
##            SBUX.close
## SBUX.close  0.2304462
mu
##            SBUX.close
## SBUX.close  0.4810674

#Calculo parámetros iniciales manera 2

x0<-tsNovoCure[1]
x <- tsNovoCure
gBm <- setModel(drift="mu*x", diffusion="sigma*x", xinit=x0)
## Warning in yuima.warn("Solution variable (lhs) not specified. Trying to use state variables."): 
## YUIMA: Solution variable (lhs) not specified. Trying to use state variables.
mod <- setYuima(model=gBm, data=setData(tsNovoCure, delta=Delta))
set.seed(123)
fit <- qmle(mod, start=list(mu=0, sigma=1),
            lower=list(mu=0.1, sigma=0.1),
            upper=list(mu=100, sigma=10))
summary(fit)
## Quasi-Maximum likelihood estimation
## 
## Call:
## qmle(yuima = mod, start = list(mu = 0, sigma = 1), lower = list(mu = 0.1, 
##     sigma = 0.1), upper = list(mu = 100, sigma = 10))
## 
## Coefficients:
##        Estimate Std. Error
## sigma 0.2330598 0.01048339
## mu    0.4813811 0.28104572
## 
## -2 log L: 702.5723
tsNovoCure[1]
## [1] 64.32
x0
## [1] 64.32

comparación

summary(fit)
## Quasi-Maximum likelihood estimation
## 
## Call:
## qmle(yuima = mod, start = list(mu = 0, sigma = 1), lower = list(mu = 0.1, 
##     sigma = 0.1), upper = list(mu = 100, sigma = 10))
## 
## Coefficients:
##        Estimate Std. Error
## sigma 0.2330598 0.01048339
## mu    0.4813811 0.28104572
## 
## -2 log L: 702.5723
coef(fit)
##     sigma        mu 
## 0.2330598 0.4813811
sigma
##            SBUX.close
## SBUX.close  0.2304462
x0
## [1] 64.32
mu
##            SBUX.close
## SBUX.close  0.4810674
gbm_vec <- function(nsim = 10000, t = 366, mu = 0.4810674, sigma = 0.2304462, S0 = 64.32, dt = 1./365) {
  # matrix of random draws - one for each day for each simulation
  epsilon <- matrix(rnorm(t*nsim), ncol = nsim, nrow = t)
  # get GBM and convert to price paths
  gbm <- exp((mu - sigma * sigma / 2) * dt + sigma * epsilon * sqrt(dt))
  gbm <- apply(rbind(rep(S0, nsim), gbm), 2, cumprod)
  return(gbm)
}
valores_simulados <- simulate(gBm, true.parameter = list(mu=mu, sigma=sigma))
## Warning in yuima.warn("'delta' (re)defined."): 
## YUIMA: 'delta' (re)defined.

#PROBAR CON FORMA 1 Y FORMA 2

#mu 1 = 0.4810674
#sigma 1 = 0.2330598
#x0 = 64.32

#mu 2 = 0.4813811
# sigma 2 = 0.2330598
#x0 = 64.32
nsim <-1000
t <-366
mu <-0.4813811
sigma <-0.2330598
S0 <-64.32
dt =1/365
nsim <-1000
t <-366
mu <-0.4813811
sigma <-0.2330598
S0 <-64.32
dt =1/365
gbm <- gbm_vec(nsim, t, mu, sigma, S0, dt)
gbm_df <- as.data.frame(gbm) %>%
  mutate(ix = 1:nrow(gbm)) %>%
  pivot_longer(-ix, names_to = 'sim', values_to = 'price')
gbm_df %>%
  ggplot(aes(x=ix, y=price, color=sim)) +
  geom_line() +
  theme(legend.position = 'none')

data.frame(price = gbm[253, ]) %>%
  ggplot(aes(x = price)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.1) +
  geom_density() + 
  ggtitle('terminal price distribution')

D <- gbm[253, ] %>%
  density()

resultado del modelo

D$x[which.max(D$y)]
## [1] 83.17347

Resultado 1 = 83.15546 Resultado 2 = 83.1013