En este documento se explica lo relacionado al pronóstico de series de tiempo. Lo cual comprende:
Elaborando la serie de tiempo en R Análisis de la serie de tiempo Pronóstico de la serie de tiempo Estimación de errores de pronóstico Para el análisis de series de tiempo en R, existen una gran variedad de librerias. En este tutorial se emplea la libreria fpp2, la cual contiene contiene:
La libreria forecast, que contiene funciones de pronóstico La libreria ggplot2, que contiene funciones para gráficos Esta libreria se instala y se carga mediante install.packages(“fpp2”)y library(fpp2) respectivamente.
Elaborando la serie de tiempo en R Cargando librerias y datos Recordemos que se puede leer datos directamente desde una API, SQL u otros medios. Para este tutorial empleo la lectura de datos desde Excel. Entonces para ello emplearemos las siguientes funciones:
La función head() nos muestra los primeros datos de nuestra tabla de datos
library(fpp2)## Warning: package 'fpp2' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──## ✔ ggplot2 3.4.2 ✔ fma 2.5
## ✔ forecast 8.22.0 ✔ expsmooth 2.3## Warning: package 'forecast' was built under R version 4.3.3## Warning: package 'fma' was built under R version 4.3.3## Warning: package 'expsmooth' was built under R version 4.3.3## AirPassengers## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1949 112 118 132 129 121 135 148 148 136 119 104 118
## 1950 115 126 141 135 125 149 170 170 158 133 114 140
## 1951 145 150 178 163 172 178 199 199 184 162 146 166
## 1952 171 180 193 181 183 218 230 242 209 191 172 194
## 1953 196 196 236 235 229 243 264 272 237 211 180 201
## 1954 204 188 235 227 234 264 302 293 259 229 203 229
## 1955 242 233 267 269 270 315 364 347 312 274 237 278
## 1956 284 277 317 313 318 374 413 405 355 306 271 306
## 1957 315 301 356 348 355 422 465 467 404 347 305 336
## 1958 340 318 362 348 363 435 491 505 404 359 310 337
## 1959 360 342 406 396 420 472 548 559 463 407 362 405
## 1960 417 391 419 461 472 535 622 606 508 461 390 432Creando el objeto de serie de tiempo La serie de tiempo se debe almacenar en el objeto ts empleando la función ts(). El empleo de esta función es de la siguiente forma : serie <- ts(data, frequency= ,start=). Donde el primer agumento es un vector de datos. Luego los dos siguientes argumentos se indican de acuerdo a la siguiente tabla.
| Frecuencia | frequency= |
start= |
|---|---|---|
| Anual | 1 | 2000 |
| Trimestral | 4 | c(2000,2) |
| Mensual | 12 | c(2000,9) |
| Diario | 7 o 365.25 | 1 o c(2000,234) |
| Semanal | 52.18 | c(2000,23) |
| Horario | 24 o 168 o 8766 | 1 |
| Cada 30 min | 48 o 336 o 17532 | 1 |
En el siguiente fragmento de código definimos el objeto de la serie de tiempo y lo almacenaremos en data_serie.
file_path <- "E:/Cambio Climatico/EXCELSTESIS/DatosTemp.csv"
dataventanastemp<- read.csv(file_path)
dataventanastemp## DATOS
## 1 24.80000
## 2 25.20000
## 3 25.60000
## 4 25.80000
## 5 24.90000
## 6 23.50000
## 7 22.70000
## 8 23.20000
## 9 25.40000
## 10 24.70000
## 11 24.70000
## 12 24.00000
## 13 24.90000
## 14 25.70000
## 15 26.20000
## 16 26.10000
## 17 24.80000
## 18 22.90000
## 19 23.00000
## 20 23.00000
## 21 25.50000
## 22 25.20000
## 23 24.60000
## 24 24.20000
## 25 25.80000
## 26 25.50000
## 27 26.20000
## 28 26.00000
## 29 26.00000
## 30 24.60000
## 31 24.00000
## 32 23.30000
## 33 25.50000
## 34 25.10000
## 35 24.90000
## 36 24.70000
## 37 25.40000
## 38 25.90000
## 39 26.40000
## 40 26.40000
## 41 25.80000
## 42 23.90000
## 43 23.20000
## 44 23.90000
## 45 25.60000
## 46 25.30000
## 47 25.10000
## 48 23.60000
## 49 26.30000
## 50 26.10000
## 51 26.40000
## 52 26.20000
## 53 25.60000
## 54 24.00000
## 55 23.40000
## 56 23.50000
## 57 26.30000
## 58 25.30000
## 59 25.00000
## 60 24.40000
## 61 26.30000
## 62 25.50000
## 63 25.90000
## 64 26.30000
## 65 25.20000
## 66 24.10000
## 67 23.60000
## 68 23.50000
## 69 25.70000
## 70 25.20000
## 71 24.30000
## 72 24.70000
## 73 25.60000
## 74 25.80000
## 75 26.50000
## 76 26.10000
## 77 25.20000
## 78 24.10000
## 79 23.50000
## 80 24.40000
## 81 25.80000
## 82 25.30000
## 83 25.30000
## 84 25.10000
## 85 25.80000
## 86 26.10000
## 87 26.10000
## 88 26.10000
## 89 25.20000
## 90 24.30000
## 91 23.70000
## 92 23.00000
## 93 24.80000
## 94 24.60000
## 95 23.90000
## 96 24.00000
## 97 24.20000
## 98 25.40000
## 99 25.90000
## 100 26.20000
## 101 25.20000
## 102 24.50000
## 103 24.30000
## 104 24.20000
## 105 25.50000
## 106 24.70000
## 107 24.20000
## 108 24.30000
## 109 25.40000
## 110 25.40000
## 111 25.90000
## 112 26.20000
## 113 25.90000
## 114 24.60000
## 115 24.30000
## 116 24.50000
## 117 26.50000
## 118 25.80000
## 119 25.20000
## 120 24.90000
## 121 26.10000
## 122 26.40000
## 123 26.70000
## 124 26.80000
## 125 26.40000
## 126 24.80000
## 127 24.20000
## 128 24.00000
## 129 24.40000
## 130 24.00000
## 131 24.70000
## 132 24.40000
## 133 25.20000
## 134 25.70000
## 135 26.50000
## 136 26.10000
## 137 25.90000
## 138 25.10000
## 139 24.60000
## 140 23.70000
## 141 26.00000
## 142 25.10000
## 143 24.20000
## 144 24.60000
## 145 24.80000
## 146 25.50000
## 147 26.30000
## 148 26.50000
## 149 26.10000
## 150 25.40000
## 151 24.10000
## 152 23.80000
## 153 25.90000
## 154 25.70000
## 155 24.70000
## 156 24.40000
## 157 25.10000
## 158 25.50000
## 159 26.20000
## 160 26.10000
## 161 25.20000
## 162 23.90000
## 163 23.00000
## 164 23.40000
## 165 26.00000
## 166 24.90000
## 167 24.90000
## 168 24.80000
## 169 25.10000
## 170 25.20000
## 171 26.50000
## 172 26.10000
## 173 26.00000
## 174 25.60000
## 175 24.80000
## 176 24.50000
## 177 26.40000
## 178 25.10000
## 179 24.70000
## 180 24.50000
## 181 25.30000
## 182 26.10000
## 183 26.40000
## 184 26.90000
## 185 26.70000
## 186 26.20000
## 187 25.50000
## 188 24.60000
## 189 26.60000
## 190 25.70000
## 191 25.70000
## 192 25.80000
## 193 26.20000
## 194 26.60000
## 195 26.90000
## 196 26.70000
## 197 26.80000
## 198 25.70000
## 199 24.80000
## 200 24.90000
## 201 26.40000
## 202 25.50000
## 203 24.80000
## 204 25.20000
## 205 25.80000
## 206 26.40000
## 207 26.30000
## 208 26.60000
## 209 26.50000
## 210 25.10000
## 211 23.90000
## 212 24.20000
## 213 25.20000
## 214 24.40000
## 215 25.10000
## 216 24.60000
## 217 25.70000
## 218 25.30000
## 219 26.00000
## 220 26.30000
## 221 25.50000
## 222 23.90000
## 223 23.50000
## 224 23.50000
## 225 25.60000
## 226 25.60000
## 227 24.40000
## 228 24.60000
## 229 25.60000
## 230 26.50000
## 231 26.40000
## 232 26.70000
## 233 26.10000
## 234 25.20000
## 235 24.10000
## 236 23.50000
## 237 25.74211
## 238 25.20000
## 239 24.50000
## 240 24.30000
## 241 26.10000
## 242 26.40000
## 243 26.26500
## 244 26.31000
## 245 25.75000
## 246 24.30000
## 247 24.00000
## 248 24.50000
## 249 25.74211
## 250 25.50000
## 251 25.20000
## 252 25.70000
## 253 25.40000
## 254 25.90000
## 255 26.30000
## 256 26.10000
## 257 25.40000
## 258 24.50000
## 259 24.00000
## 260 24.40000
## 261 25.00000
## 262 24.60000
## 263 25.20000
## 264 24.90000View(dataventanastemp)
data_serie <- ts(dataventanastemp$DATOS, frequency=12, start=2000)
View(data_serie)#Nota: Los resultados de esta guía corresponden a Data1 del archivo excel. Sugiero que tambien practiques con las otras series de tiempo como Data2 o Data3.
Análisis de la serie de tiempo Para el análisis de datos, esta libreria fpp2 nos permite emplear las siguientes funciones:
autoplot(), para graficar la serie de tiempo ggseasonplot(), para graficar la estacionalidad de una serie de tiempo ggsubseriesplot(), para graficar subseries ggAcf(), para graficar la autocorrelación decompose(), permite realizar una descomposición de estacionalidad y tendencia. Grafico de la serie de tiempo Para el gráfico simple de la serie de tiempo empleamos la funcion autoplot(). Donde solo se requiere ingresar el objeto de la serie de tiempo.
autoplot(data_serie)+
labs(title = "Serie de tiempo",
x = "Tiempo",
y = "Valor",
colour = "#00a0dc")+
theme_bw() 
Descomposición de la serie de tiempo Para la descomposición de la serie de tiempo se emplea la función decompose(). Donde se debe indicar el objeto de la serie de tiempo y el tipo de descomposición. Los tipos de descomposición que acepta esta función es additive y multiplicative.
Aditivo: Serie=T + S + I
Multiplicativo: Serie=T x S x I
Donde:
T: Tendencia
S: Estacionalidad
I: Irregular o error
# Descomposición de la serie de tiempo. Se almacena en el objeto fit
fit <- decompose(data_serie, type='additive')
#fit <- decompose(data_serie, type='multiplicative')
# Para graficar esta descomposición volvemos a emplear la funcion autoplot, pero con el objeto fit
autoplot(fit)+
labs(title = "Descomposicion de la serie de tiempo",
x = "Tiempo",
y = "Valor",
colour = "Gears")+
theme_bw()
library(highcharter)## Warning: package 'highcharter' was built under R version 4.3.3hchart(stl(data_serie, s.window='periodic'))## Warning: Deprecated function. Use the `create_axis` function.Grafico de la serie de tiempo con su tendencia El siguiente fragmento de código nos permite graficar la serie de tiempo con su tendencia. Notese que emplea el objeto fit en el cual guardamos previamente los valores de la descomposición. Nótese que se emplea la funcion trendcycle() para obtener los datos de tendencia del objeto fit.
autoplot(data_serie, series="Serie tiempo") +
autolayer(trendcycle(fit), series="Tendencia") +
labs(title = "Serie de tiempo",
x = "Tiempo",
y = "Valor"
) +
theme_bw()## Warning: Removed 12 rows containing missing values (`geom_line()`).
Grafico de estacionalidad Para realizar el gráfico de estacionalidad empleamos la función ggseasonplot. Donde el argumento es el objeto que contiene la serie de tiempo.
ggseasonplot(data_serie)
library(TSstudio)## Warning: package 'TSstudio' was built under R version 4.3.3ts_seasonal(data_serie, type = "all")ts_heatmap(data_serie)#Graficar la autocorrelacion
acf(data_serie, lag=35, col="2",main ="",
xlab="Lag",ylab="Valores")
#Graficar la autocorrelacion parcial
pacf(data_serie, lag=35, col="2",main ="",
xlab="Lag",ylab="Valores")
#Como hay confirmación de tendencia se hace la diferencia (1 lag)
plot(diff(data_serie),xlab="Tiempo",ylab="Valores")
acf(diff(data_serie), lag=35, col="2",main = "",
xlab="Lag", ylab="Valores")
pacf(diff(data_serie), lag=35,col="2",main = "",
xlab="Lag",ylab="Valores")
#Como la variancia se ve incrementa en el tiempo se debe sacar logaritmo natural de la serie
plot(diff(log(data_serie)),xlab="Tiempo",ylab="Valores")
acf(diff(log(data_serie)), lag=35, col="2",main = "",
xlab="Lag", ylab="Valores")
pacf(diff(log(data_serie)), lag=35,col="2",main = "",
xlab="Lag",ylab="Valores")
#Augmented Dickey-Fuller test para ver no estacionaridad
#Devuelve la probabilidad de ser no estacionaria
#Si es mayor que 0.05 la serie es no estacionaria (tiene tendencia)
#Si es menor que 0.05 es estacionaria (No hay tendencia)
library(aTSA)## Warning: package 'aTSA' was built under R version 4.3.2##
## Attaching package: 'aTSA'## The following object is masked from 'package:forecast':
##
## forecast## The following object is masked from 'package:graphics':
##
## identifyaTSA::adf.test(diff(log(data_serie)))## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -15.9 0.01
## [2,] 1 -13.8 0.01
## [3,] 2 -15.4 0.01
## [4,] 3 -11.7 0.01
## [5,] 4 -11.5 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -15.9 0.01
## [2,] 1 -13.7 0.01
## [3,] 2 -15.4 0.01
## [4,] 3 -11.7 0.01
## [5,] 4 -11.5 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -15.9 0.01
## [2,] 1 -13.7 0.01
## [3,] 2 -15.4 0.01
## [4,] 3 -11.7 0.01
## [5,] 4 -11.5 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01library(forecast)
#Ajuste automático de los coeficientes
ARIMAmodel <- auto.arima(log(data_serie))
ARIMAmodel## Series: log(data_serie)
## ARIMA(1,0,1)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## 0.7362 -0.2586 -0.8455 1e-04
## s.e. 0.0805 0.1196 0.0517 1e-04
##
## sigma^2 = 0.0002867: log likelihood = 664.75
## AIC=-1319.5 AICc=-1319.25 BIC=-1301.85#Ajuste automático de los coeficientes
ARIMAmodel <- auto.arima(diff(log(data_serie)))
ARIMAmodel## Series: diff(log(data_serie))
## ARIMA(0,0,2)(0,1,1)[12]
##
## Coefficients:
## ma1 ma2 sma1
## -0.4701 -0.2365 -0.8350
## s.e. 0.0611 0.0675 0.0496
##
## sigma^2 = 0.0003009: log likelihood = 655.5
## AIC=-1303.01 AICc=-1302.84 BIC=-1288.9Pronóstico Métodos simples Para el pronóstico de series de tiempo mediante métodos básicos, la libreria fpp2 nos brinda las siguientes funciones:
naive(), metodo de naive simple ses(), exponential smoothing meanf(), media movil snaive(), metodo naive considerando estacionalidad El argumento a colocar en estas funciones es la serie de tiempo y el valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar. Por ejemplo si deseamos pronosticar 12 datos, se debe indicar h=12.
Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene un ajuste con la data historica checkresiduals(), permite analizar los residuales
# elaborando el método
m1 <- snaive(data_serie, h=96)
# graficando el pronóstico
autoplot(m1)
# verificando el ajuste del método
autoplot(m1)+autolayer(fitted(m1), series="Ajuste")## Warning: Removed 12 rows containing missing values (`geom_line()`).
# verificando los residuales
checkresiduals(m1)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 213.9, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24Método holt winters Para el pronóstico de series de tiempo mediante holt winters, la libreria fpp2 nos brinda la función hw().
Los argumentos a colocar en esta funcion son:
La serie de tiempo El valor de h. Este valor de h es la cantidad de datos que deseamos pronosticar. El tipo de descomposición a usar para la estacionalidad. Los tipos de descomposición que acepta esta función es additive y multiplicative. Finalmente, para verificar el ajuste del método podemos emplear las siguientes funciones:
fitted(), obtiene pronostico con la data historica checkresiduals(), permite analizar los residuales
# elaborando el pronostico
m3 <- hw(data_serie, h=96, seasonal = 'multiplicative')
# graficando el pronóstico
autoplot(m3)
# verificando el ajuste del método
autoplot(m3)+autolayer(fitted(m3), series="Ajuste")
# verificando los residuales
checkresiduals(m3)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 54.945, df = 24, p-value = 0.0003185
##
## Model df: 0. Total lags used: 24