NFLXdata <- pdfetch_YAHOO("NFLX",from = c("2019-01-01"),to = c("2023-03-02"), interval = '1d')
Netflix <- NFLXdata[,4]
length(Netflix)
## [1] 1048
tsNetflix <- ts(Netflix, start = c(2019,1),frequency=365)
plot(tsNetflix)
Debido a que la serie de tiempo no es estacionaria, se realiza una diferenciación con logaritmos
l_Netflix<-diff(log(tsNetflix))
plot(l_Netflix)
tseries::adf.test(l_Netflix)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## Warning in tseries::adf.test(l_Netflix): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: l_Netflix
## Dickey-Fuller = -9.9501, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
Ahora se obtuvo que la serie de tiempo sí es estacionaria, ya que al realizar la prueba de Dickey-Fuller se obtuvo un p-value de 0.01 < 0.05.
Delta <- 1/365
alpha <- mean(l_Netflix)/Delta
sigma <- sqrt(var(l_Netflix)/Delta)
mu <- alpha +0.5*sigma^2
x0<-tsNetflix[1]
x <- tsNetflix
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(tsNetflix, 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.5752596 0.01257284
## mu 0.2277205 0.33965416
##
## -2 log L: 8088.01
coef(fit)
## sigma mu
## 0.5752596 0.2277205
sigma
## NFLX.close
## NFLX.close 0.5955964
mu
## NFLX.close
## NFLX.close 0.232455
Al calcular los parámetros necesarios con ambas formas, se obtuvieron valores muy cercanos. Para continuar con la predicción de los valores, se utilizará la forma 2.
gbm_vec <- function(nsim = 1048, t = 1053, mu = 0.2277205, sigma = 0.5752596, S0 = 267.66, 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)
}
nsim <- 1048
t <- 1053
mu <- 0.2277205
sigma <- 0.5752596
S0 <- 267.66
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[1053, ]) %>%
ggplot(aes(x = price)) +
geom_histogram(aes(y = ..density..), binwidth = 0.1) +
geom_density() +
ggtitle('terminal price distribution')
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
D <- gbm[1049, ] %>%
density()
D <- gbm[1050, ] %>%
density()
D <- gbm[1051, ] %>%
density()
D <- gbm[1052, ] %>%
density()
D <- gbm[1053, ] %>%
density()
D$x[which.max(D$y)]
## [1] 164.4848
Se estima un precio de 168.5434 para cada acción de Netflix para el 6 de marzo de 2023.