Load packages

library(readxl)
library(ggplot2)
library(forecast)
library(tseries)

pfbcolombia <- read_excel("C:/Users/samora/OneDrive/Estudios/Mercado de Valores/Reversión a la Media/Stocks from 2012-2016/PFBCOLOMBIA 2012-2016.xlsx", 
     sheet = "data")

Examine data

pfbcolombia$date <- as.Date(pfbcolombia$`Exchange Date`) #transformando los valores de fecha en formato de fecha.

pfbcolombia$price <- as.numeric(pfbcolombia$Close)

ggplot(pfbcolombia, aes(date, price, colour="PF Bancolombia"))+geom_line()+ggtitle("Historico de 4 años de Preferencial Bancolombia")+ylab("Precio de cierre")+xlab("Tiempo")

Ahora vamos a limpiar lo datos de la base actual.

price_ts <- ts(pfbcolombia[, c("price")])

pfbcolombia$clean_price <- tsclean(price_ts)

ggplot(pfbcolombia, aes(date, clean_price, colour="PF Bancolombia"))+geom_line()+ggtitle("Historico de 4 años de Preferencial Bancolombia")+ylab("Precio de cierre")+xlab("Tiempo")
## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.

Una vez eliminado y limpiado los outlaiers y datos vacios procederemos a usar el MA(q)

pfbcolombia$price_ma90 <- ma(pfbcolombia$clean_price, order = 90) #usando el nuevo vector de serie de tiempo, sin outliers.

pfbcolombia$price_ma275 <- ma(pfbcolombia$clean_price, order = 275) #número de días de negociación en un año.

ggplot()+geom_line(data = pfbcolombia, aes(x= date, y=clean_price, colour="Precio de cierre"))+geom_line(data=pfbcolombia, aes(x=date, y=price_ma90, colour="Promedio Móvil Trimestral"))+geom_line(data = pfbcolombia, aes(x=date, y=price_ma275, colour="Promedio Móvil Anual"))+ggtitle("Preferencial Bancolombia")+ylab("Precio de Cierre")+xlab("Tiempo")
## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
## Warning: Removed 90 rows containing missing values (geom_path).
## Warning: Removed 274 rows containing missing values (geom_path).

Ahora vamos a descomponer la serie con el fin de encontrar el componente estacional, el ciclo y el error.

price_ma <- ts(na.omit(pfbcolombia$price_ma90), frequency = 30)

decomp <- stl(price_ma, s.window = "periodic")

deseasonal_price <- seasadj(decomp)

plot(decomp)

Estacionalidad

Inicaremos con el Test aumentado de Dickey - Fuller de Estacionalidad, con este podemos rechazar la hipotesis alternativa de estacionalidad de la serie de precios de bancolombia preferencial.

adf.test(price_ma, alternative = "stationary")
## Warning in adf.test(price_ma, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  price_ma
## Dickey-Fuller = -3.9847, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

Autocorrelación y Elección en el orden del Modelo

library(forecast)

Acf(price_ma, main="Autocorrelograma")

Pacf(price_ma, main="Autocorrelograma Parcial")

Ahora revisaremos las primeras diferencias y miraremos si es necesaria aplicarse.

price_d1 <- diff(deseasonal_price, differences=1)

plot(price_d1)

adf.test(price_d1, alternative = "stationary")
## Warning in adf.test(price_d1, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  price_d1
## Dickey-Fuller = -4.0399, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary

Ya bien, teniendo la seríe con primeras diferencias podemos mejorar la observación de retrasos o “lag” dentro de la Autocorrelación. Esto nos sugiere así mismo el orden a implementar en los componentes de Promedio Móvil (MA) y de Autocorrelación (AR).

Acf(price_d1, main="Autocorrelograma para la Serie de Primeras Diferencias")

Pacf(price_d1, main="Autocorrelograma Parcial para la Serie de Primeras Diferencias")

Con Primeras diferencias podemos ver que las innovaciones aún sobrepasan los intervalos de confianza, razón por la cual realizaremos el proceso con segundas diferencias y ver si mejorá el proceso para observar los retrasos en la serie de tiempo.

price_d2 <- diff(deseasonal_price, differences=2)

plot(price_d2)

adf.test(price_d2, alternative = "stationary")
## Warning in adf.test(price_d2, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  price_d2
## Dickey-Fuller = -10.754, Lag order = 11, p-value = 0.01
## alternative hypothesis: stationary
Acf(price_d2, main="Autocorrelograma para la Serie de segundas Diferencias")

Pacf(price_d2, main="Autocorrelograma Parcial para la Serie de Segundas Diferencias")

Despues de realizar el proceso con Segundas diferencias podemos encontrar entonces mayor fuerza de ruido blanco y los correlogramas encuentran un mejor nivel de observacion de los retrasos en la serie del activo.

Ajustando un Modelo ARMA

Con el paquete de “Forecast” en R es posible realizar la función del modelo ARIMA de forma mas optima, ya que la funcion “auto.arima()” se busca a través de conbinaciones el orden de los parametros que optimicen y ajusten el criterio.

A pesar de que la función auto.arima() es muy eficiente es necesario realizar los pasos anteriores con el fin de que se pueda entender la serie e interpretar los resultados del modelo, por defecto el modelo nos brinda como máximo un orden de 5 en los retraso observados dentro de su proceso de optimización.

Cabe resaltar que existe un número considerable de la calidad de ajuste a través de diferentes modelos. Dos de los más usados son el Akaike Information Criteria (AIC) y el Bayesian Information Criteria (BIC).

Al elegir el modelo podemos compara dichos criterios y elegir uno que minimice el AIC y el BIC.

auto.arima(deseasonal_price, seasonal = FALSE)
## Series: deseasonal_price 
## ARIMA(3,1,1) 
## 
## Coefficients:
##          ar1     ar2      ar3     ma1
##       0.8985  0.1714  -0.0937  0.8585
## s.e.  0.0332  0.0504   0.0314  0.0184
## 
## sigma^2 estimated as 11.51:  log likelihood=-3594.27
## AIC=7198.54   AICc=7198.58   BIC=7224.67

Usando ARIMA podemos denotar que el modelo ajustado puede ser escrito como:

*** Ydt = 0.8985Yt-1 + 0.1714Yt-2 -0.0937Yt-3 + 0.8585et-1+E ***

El coeficiente AR (1) = 0.8985 nos dice que el siguiente valor de la serie se toma como un valor anterior atenuado por un factor de 0,85 y depende del retraso de error anterior.

Evaluar e Iterar

Entonces, ahora hemos ajustado un modelo que puede producir un pronóstico, pero ¿tiene sentido? ¿Podemos confiar en este modelo? Podemos comenzar examinando los gráficos ACF y PACF para los residuos del modelo. Si los parámetros de orden del modelo y la estructura se especifican correctamente, no esperaríamos que existieran autocorrelaciones significativas.

fit <- auto.arima(deseasonal_price, seasonal = FALSE)
tsdisplay(residuals(fit),  lag.max = 30, main = "(3,1,1) Residuos del Modelo")

Hay un patrón claro presente en ACF / PACF y las parcelas residuales del modelo se repiten en el desfase 28. Esto sugiere que nuestro modelo puede estar mejor con una especificación diferente, como AR(p) = 28 o MA(q) = 28.

Podemos repetir el proceso de ajuste que permite el componente MA (28) y examinar los diagramas de diagnóstico nuevamente. Esta vez, no hay autocorrelaciones significativas presentes. Si el modelo no se especifica correctamente, eso generalmente se reflejará en residuos en forma de tendencias, enredos o cualquier otro patrón no capturado por el modelo. Idealmente, los residuos deberían verse como ruido blanco, lo que significa que normalmente se distribuyen normal.

fit2 <- arima(deseasonal_price, order = c(3,1,29))
## Warning in log(s2): Se han producido NaNs
## Warning in arima(deseasonal_price, order = c(3, 1, 29)): possible
## convergence problem: optim gave code = 1
fit2
## 
## Call:
## arima(x = deseasonal_price, order = c(3, 1, 29))
## 
## Coefficients:
##          ar1     ar2     ar3     ma1     ma2    ma3     ma4     ma5
##       0.2617  0.2067  0.4654  1.4928  1.0844  0.465  0.0774  0.0694
## s.e.  0.3273  0.0925  0.2665  0.3158  0.5150  0.225  0.0774  0.0734
##          ma6     ma7     ma8     ma9    ma10    ma11    ma12    ma13
##       0.0389  0.0175  0.0449  0.0545  0.0621  0.0725  0.0764  0.0521
## s.e.  0.0736  0.0735  0.0730  0.0718  0.0732  0.0723  0.0739  0.0747
##         ma14     ma15     ma16    ma17    ma18     ma19     ma20    ma21
##       0.0098  -0.0108  -0.0136  0.0243  0.0066  -0.0268  -0.0397  0.0115
## s.e.  0.0718   0.0730   0.0755  0.0750  0.0723   0.0721   0.0752  0.0724
##         ma22    ma23    ma24    ma25    ma26    ma27    ma28    ma29
##       0.0221  0.0656  0.0386  0.0367  0.0691  0.1166  0.1507  0.0538
## s.e.  0.0724  0.0706  0.0675  0.0646  0.0611  0.0606  0.0548  0.0334
## 
## sigma^2 estimated as 10.72:  log likelihood = -3582.02,  aic = 7230.05
tsdisplay(residuals(fit2), main = "Residuos del Modelo Estacional")

Revisar esta parte!!!!

La previsión o prónostico utilizando un modelo ajustado es sencilla en R. Podemos especificar el horizonte de pronóstico h periodos por delante para que se realicen las predicciones, y usar el modelo ajustado para generar esas predicciones:

fcast <- forecast(fit2, h=30)
plot(fcast)

La línea azul clara de arriba muestra el ajuste proporcionado por el modelo, pero ¿qué pasaría si quisiéramos tener una idea de cómo funcionará el modelo en el futuro? Un método es reservar una parte de nuestros datos como un conjunto de “Hold-out”, ajustar el modelo y luego comparar el pronóstico con los valores reales observados:

hold <- window(ts(deseasonal_price), start=1300)
fit_no_holdout <- arima(ts(deseasonal_price[-c(1300:1375)]), order =c(3,1,29))
## Warning in log(s2): Se han producido NaNs

## Warning in log(s2): Se han producido NaNs

## Warning in log(s2): Se han producido NaNs
## Warning in arima(ts(deseasonal_price[-c(1300:1375)]), order = c(3, 1, 29)):
## possible convergence problem: optim gave code = 1
fcast_no_holdout <- forecast(fit_no_holdout, h=75)
plot(fcast_no_holdout, main = " ")
lines(ts(deseasonal_price))

Podemos ver que la linea azul, la cual representa el pronostico, es es una linea recta,lo cual parece inprobable dado el comportamiento anerior. Es necesario recordar que este modelo no se asume estacionalidad.

Vamos ahora a asignar un modelo con estacionalidad, con el fin de ver como mejora la predicción del precio en la acción.

fit_w_seasonality <- auto.arima(deseasonal_price, seasonal = TRUE)
fit_w_seasonality
## Series: deseasonal_price 
## ARIMA(3,1,1)(0,0,2)[30] 
## 
## Coefficients:
##          ar1     ar2      ar3     ma1    sma1    sma2
##       1.0102  0.0056  -0.0508  0.8956  0.4931  0.4927
## s.e.  0.0307  0.0471   0.0298  0.0132  0.0624  0.0615
## 
## sigma^2 estimated as 10.39:  log likelihood=-3528.91
## AIC=7071.83   AICc=7071.91   BIC=7108.41
seas_fcast <- forecast(fit_w_seasonality, h=75)
plot(seas_fcast)

Ahora revisemos los residuos:

fit3 <- arima(deseasonal_price, order = c(3,1,29), seasonal = list(order=c(0,0,2), period=30))
## Warning in arima(deseasonal_price, order = c(3, 1, 29), seasonal =
## list(order = c(0, : possible convergence problem: optim gave code = 1
fit3
## 
## Call:
## arima(x = deseasonal_price, order = c(3, 1, 29), seasonal = list(order = c(0, 
##     0, 2), period = 30))
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): Se han producido NaNs
##          ar1      ar2     ar3     ma1     ma2     ma3     ma4    ma5
##       1.2442  -0.8705  0.2956  0.4575  0.3274  0.8834  0.3689  0.596
## s.e.     NaN      NaN  0.0149     NaN     NaN  0.0315     NaN    NaN
##          ma6     ma7     ma8     ma9    ma10    ma11    ma12    ma13
##       0.6311  0.4134  0.6064  0.5445  0.4301  0.6689  0.4519  0.4786
## s.e.  0.0348     NaN  0.0529     NaN     NaN  0.0637     NaN     NaN
##         ma14    ma15   ma16    ma17    ma18    ma19    ma20    ma21
##       0.6884  0.3894  0.579  0.6428  0.3779  0.6346  0.5636  0.4463
## s.e.  0.0589     NaN    NaN  0.0360     NaN  0.0552     NaN     NaN
##         ma22    ma23    ma24    ma25    ma26    ma27    ma28    ma29
##       0.6678  0.4745  0.5239  0.6505  0.4066  0.6316  0.6763  0.2012
## s.e.  0.0782     NaN     NaN  0.0761     NaN     NaN     NaN     NaN
##         sma1    sma2
##       0.9241  0.8986
## s.e.  0.0212  0.0270
## 
## sigma^2 estimated as 7.169:  log likelihood = -3361.06,  aic = 6792.11
tsdisplay(residuals(fit3), main = "Residuos del Modelo Estacional")

Realicemos entonces el prnostico con este nuevo ajuste y observemos que los el modelo tiene un desfase de 25. Lo que nos sugiere un cambio en el ajuste, sin embargo veremos el pronóstico y ver su mejora en comparación con el ajuste anterior:

hold <- window(ts(deseasonal_price), start=1300)
fit_no_holdout3 <- arima(ts(deseasonal_price[-c(1300:1375)]), order = c(3,1,29), seasonal = list(order=c(0,0,2), period=30))
## Warning in log(s2): Se han producido NaNs

## Warning in log(s2): Se han producido NaNs
## Warning in arima(ts(deseasonal_price[-c(1300:1375)]), order = c(3, 1,
## 29), : possible convergence problem: optim gave code = 1
fcast_no_holdout3 <- forecast(fit_no_holdout3, h=75)
plot(fcast_no_holdout3, main = " ")
lines(ts(deseasonal_price))

Ahora mirando la predicción únicamente con los últimos 30 días.

hold <- window(ts(deseasonal_price), start=1345)
fit_no_holdout3 <- arima(ts(deseasonal_price[-c(1345:1375)]), order = c(3,1,29), seasonal = list(order=c(0,0,2), period=30))
## Warning in log(s2): Se han producido NaNs
fcast_no_holdout3 <- forecast(fit_no_holdout3, h=30)
plot(fcast_no_holdout3, main = " ")
lines(ts(deseasonal_price))