Deseamos pronosticar el desempeño de la marca Apple. Para eso instalaremos las librerías.

#library(pdfetch)
#library(tidyverse)
#library(yuima)

Cargamos los datos de Apple, estos se obtuvieron de Yahoo!finanzas. Se cargaran en total los datos de 4 años (2019-2023).

AAPLdata <- pdfetch_YAHOO("AAPL",from = c("2019-01-01"),to = c("2023-02-18"), interval = '1d')

#Obtenemos la columna que es de nuestro interés
Apple <- AAPLdata[,4]
length(Apple)
## [1] 1041
#Convertimos en una serie de tiempo
tsApple <- ts(Apple, start = c(2019,1),frequency=365)

Gráfico de nuestra serie de datos

Ver si es estacionaria

plot(tsApple)

tseries::adf.test(tsApple)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
##  Augmented Dickey-Fuller Test
## 
## data:  tsApple
## Dickey-Fuller = -1.819, Lag order = 10, p-value = 0.6549
## alternative hypothesis: stationary

Como se puede observar en el gráfico, la serie de datos no es estacionaria.Ya que esta presenta variaciones como lo son crecimientos y dismiuciones en el gráfico. También el p-value es mayor a 0.946. Para lograr que esta sea estacionaria debemos que realizar una diferenciación.

Primera diferenciación

l_Apple<-diff(log(tsApple))
plot(l_Apple)

tseries::adf.test(l_Apple)
## Warning in tseries::adf.test(l_Apple): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  l_Apple
## Dickey-Fuller = -9.4941, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary

Luego de realizar la primera diferenciación pudimos obtener un p-value menor a 0.05.

Parámetros iniciales

Manera 1

Delta <- 1/365
alpha <- mean(l_Apple)/Delta
sigma <- sqrt(var(l_Apple)/Delta)
mu <- alpha +0.5*sigma^2
x0<-tsApple[1]

Manera 2

x <- tsApple
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(tsApple, 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.4142429 0.009114667
## mu    0.5601307 0.245405804
## 
## -2 log L: 4592.968

Comparamos los valores obtenidos de la manera 1 y 2

coef(fit)
##     sigma        mu 
## 0.4142429 0.5601307
sigma
##            AAPL.close
## AAPL.close  0.4136048
mu
##            AAPL.close
## AAPL.close  0.5599285
gbm_vec <- function(nsim = 10000, t = 25, mu = 0, sigma = 0.1, S0 = 100, 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)
}

Gráfico de las acciones

nsim <- 1000
t <- 1046
mu <- 0.5186177
sigma <- 0.4164425
S0 <- x0
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')

Como se puede ver el el gráfico el precio de la acción de Apple va en aumento. Esto quiere decir que le va muy bien.

Gráfico de la distribución del precio

data.frame(price = gbm[253, ]) %>%
  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[253, ] %>%
  density()

#resultado del modelo
D$x[which.max(D$y)]  
## [1] 49.05642
NVCRdata2 <- pdfetch_YAHOO("AAPL",from = c("2019-01-01"),to = c("2023-02-18"), interval = '1d')
Apple2 <- NVCRdata2[,4]
View(Apple2)