SESIÓN 5. Pronósticos: Invertir diferencias regulares en R

Análisis de Series Temporales

Activar librerías

#install.packages("devtools")
#devtools::install_github("r-lib/conflicted")
library(conflicted)
library(readxl)
library(dplyr)
library(stats)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ ggplot2   3.5.1     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1

library(tseries)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

library(urca)
library(forecast)

Dataset de comodity_price

library(readr)
comodity_price <- read_csv("comodity_price.csv")

Rows: 755 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl  (6): Open, High, Low, Close, Adj Close, Volume
date (1): Date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(comodity_price)
dim(comodity_price)

[1] 755   7

gold_price=ts(comodity_price$Open)
library(TSstudio)
ts_plot(gold_price)

tail(gold_price)

Time Series:
Start = 750 
End = 755 
Frequency = 1 
[1] 16.55 16.51 16.46 16.24 16.56 16.61

• Definir si cumple o no con estacionariedad

adf.test(gold_price)

    Augmented Dickey-Fuller Test

data:  gold_price
Dickey-Fuller = -2.7386, Lag order = 9, p-value = 0.2657
alternative hypothesis: stationary

• Como el p\_valor es mayor que 0.05, entonces la serie es no estacionaria

tsdisplay(gold_price)

Diferenciamos la serie

• La diferencia regular, ayuda a que la serie sea estacionaria

gold_price_diff <- diff(gold_price, lag = 1)
tm <- cbind(gold_price, gold_price_diff)
head(tm)

Time Series:
Start = 1 
End = 6 
Frequency = 1 
  gold_price gold_price_diff
1      27.10              NA
2      26.97       -0.130001
3      26.61       -0.359998
4      27.13        0.519998
5      29.62        2.490002
6      31.20        1.580000

plot(gold_price_diff)

• Definir si cumple o no con estacionariedad

adf.test(gold_price_diff)

Warning in adf.test(gold_price_diff): p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  gold_price_diff
Dickey-Fuller = -8.5322, Lag order = 9, p-value = 0.01
alternative hypothesis: stationary

tsdisplay(gold_price_diff)

• Como el p\_valor es menor que 0.05, entonces la serie es estacionaria

Modelo propuesto ARIMA(p,d,q)=ARIMA(1,1,1) o identificar en \pm 1 de AR o MA

Estimación del modelo

\underline{NOTA}. Ya es estacionaria (con la diferencia regular realizada), entonces fijamos la integración (d) en cero; es decir: “max.d=0”

myarima=auto.arima(gold_price_diff, trace = T, 
                   stepwise = F, 
                   approximation = F,max.d=0)

 ARIMA(0,0,0) with zero mean     : 950.5442
 ARIMA(0,0,0) with non-zero mean : 951.8461
 ARIMA(0,0,1) with zero mean     : 952.5455
 ARIMA(0,0,1) with non-zero mean : 953.8465
 ARIMA(0,0,2) with zero mean     : 950.1335
 ARIMA(0,0,2) with non-zero mean : 951.3222
 ARIMA(0,0,3) with zero mean     : 949.8159
 ARIMA(0,0,3) with non-zero mean : 951.111
 ARIMA(0,0,4) with zero mean     : 951.7984
 ARIMA(0,0,4) with non-zero mean : 953.0801
 ARIMA(0,0,5) with zero mean     : 952.7529
 ARIMA(0,0,5) with non-zero mean : 953.947
 ARIMA(1,0,0) with zero mean     : 952.5469
 ARIMA(1,0,0) with non-zero mean : 953.8489
 ARIMA(1,0,1) with zero mean     : 953.8358
 ARIMA(1,0,1) with non-zero mean : 955.0082
 ARIMA(1,0,2) with zero mean     : 950.4508
 ARIMA(1,0,2) with non-zero mean : 951.6769
 ARIMA(1,0,3) with zero mean     : 951.8228
 ARIMA(1,0,3) with non-zero mean : 953.1157
 ARIMA(1,0,4) with zero mean     : 952.9792
 ARIMA(1,0,4) with non-zero mean : 954.1276
 ARIMA(2,0,0) with zero mean     : 950.2278
 ARIMA(2,0,0) with non-zero mean : 951.4184
 ARIMA(2,0,1) with zero mean     : 950.1936
 ARIMA(2,0,1) with non-zero mean : 951.4224
 ARIMA(2,0,2) with zero mean     : Inf
 ARIMA(2,0,2) with non-zero mean : Inf
 ARIMA(2,0,3) with zero mean     : Inf
 ARIMA(2,0,3) with non-zero mean : Inf
 ARIMA(3,0,0) with zero mean     : 949.7236
 ARIMA(3,0,0) with non-zero mean : 951.0113
 ARIMA(3,0,1) with zero mean     : 951.7217
 ARIMA(3,0,1) with non-zero mean : 953.0067
 ARIMA(3,0,2) with zero mean     : Inf
 ARIMA(3,0,2) with non-zero mean : Inf
 ARIMA(4,0,0) with zero mean     : 951.6992
 ARIMA(4,0,0) with non-zero mean : 952.9754
 ARIMA(4,0,1) with zero mean     : 952.9343
 ARIMA(4,0,1) with non-zero mean : 954.0738
 ARIMA(5,0,0) with zero mean     : 953.1968
 ARIMA(5,0,0) with non-zero mean : 954.4156



 Best model: ARIMA(3,0,0) with zero mean     

• Se identifica como mejor modelo el ARIMA(3,1,0)

Ejecutar predicciones del modelo

arimafore <- forecast(myarima, h = 20)
plot(arimafore)

arimafore$mean

Time Series:
Start = 756 
End = 775 
Frequency = 1 
 [1] -3.685182e-02  1.469286e-02  5.687281e-03 -3.235833e-03  4.171231e-04
 [6]  5.736213e-04 -2.181138e-04 -1.940108e-05  4.962591e-05 -1.110319e-05
[11] -4.880035e-06  3.703886e-06 -2.701699e-07 -5.621844e-07  2.340831e-07
[16]  2.704686e-08 -5.016464e-08  1.144843e-08  5.363077e-09 -3.760708e-09

Realizamos la inversión

• Implica anular la diferenciación regular, para establecer el proceso con la serie original y con ello, garantizar que los pronósticos no tengan sesgo.

last_gp=gold_price[775]

pred=arimafore$mean
de_diff_gold_fore=cumsum(pred)+ last_gp

Convertir a ts

• La idea es, convertirlo a un objeto ts que empiece en el índice después del último de los precios

ddiff_gfore_ts=ts(de_diff_gold_fore,start=776)
ts.plot(gold_price, ddiff_gfore_ts, lty = c(1,3), col=c(5,2))

#colors()

#plot(ddiff_gfore_ts)

Análisis de Red Neuronal

Activar las librerías

library(forecast)
library(nnet)
library(ggplot2)
mod_nnet=nnetar(gold_price)
mod_nnet

Series: gold_price 
Model:  NNAR(7,4) 
Call:   nnetar(y = gold_price)

Average of 20 networks, each of which is
a 7-4-1 network with 37 weights
options were - linear output units 

sigma^2 estimated as 0.1727

summary(mod_nnet)

          Length Class        Mode     
x         755    ts           numeric  
m           1    -none-       numeric  
p           1    -none-       numeric  
P           1    -none-       numeric  
scalex      2    -none-       list     
size        1    -none-       numeric  
subset    755    -none-       numeric  
model      20    nnetarmodels list     
nnetargs    0    -none-       list     
fitted    755    ts           numeric  
residuals 755    ts           numeric  
lags        7    -none-       numeric  
series      1    -none-       character
method      1    -none-       character
call        2    -none-       call     

Ejecutar predicciones del modelo (red neuronal)

nnetforecast <- forecast(mod_nnet, h = 100, PI = F)
nnetforecast$mean

Time Series:
Start = 756 
End = 855 
Frequency = 1 
  [1] 16.58411 16.58567 16.56427 16.54272 16.53361 16.51873 16.50241 16.48699
  [9] 16.47079 16.45460 16.43870 16.42275 16.40682 16.39096 16.37513 16.35937
 [17] 16.34369 16.32809 16.31258 16.29716 16.28186 16.26667 16.25160 16.23666
 [25] 16.22187 16.20723 16.19274 16.17841 16.16426 16.15029 16.13650 16.12290
 [33] 16.10950 16.09630 16.08332 16.07055 16.05800 16.04567 16.03357 16.02170
 [41] 16.01007 15.99867 15.98751 15.97660 15.96593 15.95550 15.94532 15.93538
 [49] 15.92569 15.91625 15.90705 15.89809 15.88937 15.88090 15.87266 15.86466
 [57] 15.85689 15.84935 15.84204 15.83495 15.82808 15.82143 15.81499 15.80876
 [65] 15.80274 15.79691 15.79129 15.78585 15.78060 15.77554 15.77065 15.76594
 [73] 15.76139 15.75701 15.75280 15.74873 15.74482 15.74105 15.73743 15.73394
 [81] 15.73059 15.72737 15.72427 15.72129 15.71843 15.71568 15.71304 15.71051
 [89] 15.70808 15.70574 15.70351 15.70136 15.69930 15.69732 15.69542 15.69361
 [97] 15.69186 15.69019 15.68859 15.68706

autoplot(nnetforecast)