Modelo SARIMA para la series Renta
Autores:
Jessica Lopez Mejia
Marlon Orrego Corzo
Nicolas Romero
Solano
library(urca)
library(forecast)
library(tseries)
library(lmtest)
library(uroot)
library(readxl)
library(fUnitRoots)
library(tsoutliers)
library(sarima)
library(aTSA)
library(fpp)
require("PolynomF")setwd("C:/Users/Embag/OneDrive/Documents")Renta <- readxl::read_excel("C:/Users/Embag/OneDrive/Documents/datos-1.xlsx", col_names = TRUE)
Renta <- stats::ts(Renta$I_Renta, start = c(2000,01,01), frequency = 12)
plot(Renta)
Como lo vimos en la parte descriptiva, a nuestra serie sobre el Impuesto de renta se le debe estabilizar la varianza dado que en la prueba de Box-Cox obtuvimos un valor de \(\lambda=0\), lo que nos lleva a aplicarle una transformación logarítmica a nuestra serie.
forecast::BoxCox.lambda(Renta, method ="loglik", lower = -1, upper = 3)## [1] 0MASS::boxcox(lm(Renta ~ 1),seq(-2, 3, length = 50))
plot(forecast::BoxCox(Renta,lambda = 0))
# Aplicar log para estabilizar la varianza
lRenta <- log(Renta)
# Realizamos de nuevo la prueba de Box-Cox
forecast::BoxCox.lambda(lRenta, method ="loglik", lower = -5, upper = 3)## [1] 0.9MASS::boxcox(lm(lRenta ~ 1),seq(-5, 2, length = 50)) 
Ahora, el uno se encuentra en el intervalo y, por lo tanto, se ha estabilizado la varianza de la serie.
Estudiando la tendencia de la Serie Renta
El siguiente paso es detectar la presencia o no de raices unitarias. Para ello procedemos a realizar la prueba de raices unitarias de Dickey-Fuller cuyo sitema de hipótesis se define de la siguiente forma:
\[ \begin{cases} H_0: & \text{Hay presencia de raíz unitaria} \\ H_1: & \text{No hay presencia de raíz unitaria } \end{cases} \]
ar(lRenta) # Número de retardos que se pueden tener en cuenta en la prueba##
## Call:
## ar(x = lRenta)
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.3117 0.5701 0.1098 -0.0249 -0.0681 -0.0995 -0.0099 0.0816
## 9 10 11 12 13 14
## 0.0370 0.0581 -0.1125 0.5920 -0.2561 -0.2283
##
## Order selected 14 sigma^2 estimated as 0.08034fUnitRoots::adfTest(lRenta,lags = 14,type='nc')## Warning in fUnitRoots::adfTest(lRenta, lags = 14, type = "nc"): p-value greater
## than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 14
## STATISTIC:
## Dickey-Fuller: 4.6335
## P VALUE:
## 0.99
##
## Description:
## Wed Sep 25 12:49:36 2024 by user: EmbagaTSA::adf.test(lRenta,nlag = 15)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -0.148 0.601
## [2,] 1 0.641 0.828
## [3,] 2 0.947 0.907
## [4,] 3 0.889 0.900
## [5,] 4 0.853 0.889
## [6,] 5 0.955 0.908
## [7,] 6 1.120 0.930
## [8,] 7 2.014 0.990
## [9,] 8 2.359 0.990
## [10,] 9 4.205 0.990
## [11,] 10 3.496 0.990
## [12,] 11 8.134 0.990
## [13,] 12 7.202 0.990
## [14,] 13 5.256 0.990
## [15,] 14 4.633 0.990
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -7.4268 0.010
## [2,] 1 -2.2452 0.231
## [3,] 2 -1.9032 0.366
## [4,] 3 -1.8078 0.403
## [5,] 4 -1.8198 0.398
## [6,] 5 -1.6613 0.461
## [7,] 6 -1.6150 0.479
## [8,] 7 -0.9701 0.708
## [9,] 8 -0.8686 0.743
## [10,] 9 -0.4416 0.893
## [11,] 10 -0.6038 0.836
## [12,] 11 -0.2454 0.925
## [13,] 12 -0.1821 0.934
## [14,] 13 -0.0519 0.951
## [15,] 14 0.0253 0.958
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -24.78 0.0100
## [2,] 1 -7.97 0.0100
## [3,] 2 -6.64 0.0100
## [4,] 3 -7.14 0.0100
## [5,] 4 -8.00 0.0100
## [6,] 5 -8.15 0.0100
## [7,] 6 -8.45 0.0100
## [8,] 7 -5.41 0.0100
## [9,] 8 -4.95 0.0100
## [10,] 9 -2.83 0.2265
## [11,] 10 -3.49 0.0434
## [12,] 11 -1.48 0.7956
## [13,] 12 -1.42 0.8207
## [14,] 13 -1.58 0.7548
## [15,] 14 -1.62 0.7366
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01Ahora, hacemos diferenciación ordinaria para quitarle la tendencia a la serie y volvemos a aplicar la prueba de Dickey-Fuller para saber si se necesitan hacerle más diferenciaciones a la serie.
dlRenta <- diff(lRenta)
plot(dlRenta)
fUnitRoots::adfTest(dlRenta,lags = 14,type='nc')## Warning in fUnitRoots::adfTest(dlRenta, lags = 14, type = "nc"): p-value
## smaller than printed p-value##
## Title:
## Augmented Dickey-Fuller Test
##
## Test Results:
## PARAMETER:
## Lag Order: 14
## STATISTIC:
## Dickey-Fuller: -3.8923
## P VALUE:
## 0.01
##
## Description:
## Wed Sep 25 12:49:37 2024 by user: EmbagaTSA::adf.test(dlRenta,nlag = 15)## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -57.82 0.01
## [2,] 1 -20.88 0.01
## [3,] 2 -13.01 0.01
## [4,] 3 -10.04 0.01
## [5,] 4 -9.28 0.01
## [6,] 5 -8.72 0.01
## [7,] 6 -12.51 0.01
## [8,] 7 -11.17 0.01
## [9,] 8 -15.41 0.01
## [10,] 9 -9.05 0.01
## [11,] 10 -16.30 0.01
## [12,] 11 -9.38 0.01
## [13,] 12 -5.77 0.01
## [14,] 13 -4.66 0.01
## [15,] 14 -3.89 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -57.78 0.01
## [2,] 1 -20.91 0.01
## [3,] 2 -13.05 0.01
## [4,] 3 -10.08 0.01
## [5,] 4 -9.34 0.01
## [6,] 5 -8.81 0.01
## [7,] 6 -12.75 0.01
## [8,] 7 -11.52 0.01
## [9,] 8 -16.44 0.01
## [10,] 9 -9.89 0.01
## [11,] 10 -19.95 0.01
## [12,] 11 -12.55 0.01
## [13,] 12 -8.03 0.01
## [14,] 13 -6.71 0.01
## [15,] 14 -5.69 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -57.68 0.01
## [2,] 1 -20.87 0.01
## [3,] 2 -13.03 0.01
## [4,] 3 -10.06 0.01
## [5,] 4 -9.32 0.01
## [6,] 5 -8.79 0.01
## [7,] 6 -12.72 0.01
## [8,] 7 -11.50 0.01
## [9,] 8 -16.41 0.01
## [10,] 9 -9.87 0.01
## [11,] 10 -19.91 0.01
## [12,] 11 -12.53 0.01
## [13,] 12 -8.01 0.01
## [14,] 13 -6.70 0.01
## [15,] 14 -5.68 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01Con base en lo obtenido anteriormente podemos concluir que nuestra serie únicamente requiere de una diferenciación para remover la tendencia y entonces sería razonable considerar un modelo ARIMA de orden \(d=1\) para modelar nuestra serie de tiempo.
Estudiando la estacionaidad de la Serie Renta
Como lo vimos en el análisis descriptivo de nuestra serie, existe una componente estacional de periodo 12, es decir, que cada 12 meses el valor medio del impuesto de Renta se comporta de manera parecida. Podemos observar este detalle en el siguiente gráfico:
ggseasonplot(dlRenta)
Otro gráfico que nos puede mostrar si el valor medio del impuesto de renta varía a lo largo de los meses es el gráfico del valor del impuesto de renta agrupado por mes.
monthplot(dlRenta)
Gráficamente podemos deducir que nuestra serie presenta una estacionalidad anual. Sim embargo, podemos verificar este hallazgo realizando la prueba de raices estacionales .
nsdiffs(dlRenta)## [1] 1En base al resultado de la prueba procedemos a realizarle una diferenciación estacional para remover la componente estacional y luego realizamos nuevamente la prueba nsdiffs para comprobar que efectivamente hemos removido la componente estacional de nuestra serie.
DdlRenta=diff(dlRenta,lag=12) # Periodo 12
nsdiffs(DdlRenta)## [1] 0Ahora, observemos gráficamente la serie para comprobar que hemos removido la componente estacional anual de nuestra serie de tiempo.
monthplot(DdlRenta)
Teniendo en cuenta todo el análisis hecho hasta ahora podemos concluir que nuestra serie podría ser modelada con un modelo SARIMA de órden \(D=1\) y \(d=1\). Ahora, para conocer los órdenes \(p,q\) y \(P,Q\) nos basaremos en la función de autocorrelación simple ACF y la función de autocorrelación parcial PACF.
Identificando los órdenes p,q y P,Q del modelo SARIMA
Analicemos los gráficos de autocorrelación simple y parcial de nuestra serie de tiempo para determinar los órdenes de nuestro modelo SARIMA
acf(DdlRenta)
acf(DdlRenta,lag.max = 50, ci.type='ma') #max q=2, max Q=2.
pacf(DdlRenta,lag.max = 50) #max p=5, max P=0.
Con lo anterior podemos postular un modelo SARIMA(5,1,2)x(0,1,2). Sin embargo, veamos qué modelo obtenemos cuando lo buscamos de forma automática.
###Arima Automático
modelo.automatico=auto.arima(DdlRenta,d=1,D=1,max.p=11,max.q=2,start.p=0, start.q=0,seasonal=TRUE,max.order=12,stationary=TRUE,ic="aicc",stepwise=FALSE,allowmean = TRUE)
modelo.automatico## Series: DdlRenta
## ARIMA(1,0,2)(0,0,2)[12] with zero mean
##
## Coefficients:
## ar1 ma1 ma2 sma1 sma2
## -0.6169 -0.4306 -0.3132 -0.4264 -0.3060
## s.e. 0.1483 0.1674 0.1553 0.0610 0.0614
##
## sigma^2 = 0.03126: log likelihood = 84.01
## AIC=-156.03 AICc=-155.71 BIC=-134.33Bajo el modelo automático mantenemos el valor de los órdenes \(q\), \(Q\) y \(P\) pero nos aconseja un valor de \(p=1\).
Ajustando el modelo SARIMA
Ajustemos un modelo SARIMA con base en los hallazgos encontrados en la función de autocorrelación simple y autocorrelación parcial y luego analicemos cuales parámetros resultan o no significativos al \(5\%\).
fit1<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.135223 0.089828 -12.6378 < 2.2e-16 ***
## ar2 0.040693 0.152103 0.2675 0.78906
## ar3 0.237035 0.143707 1.6494 0.09906 .
## ar4 0.027225 0.124854 0.2181 0.82739
## ar5 -0.023839 0.071705 -0.3325 0.73955
## ma1 0.097781 0.067918 1.4397 0.14995
## ma2 -0.845705 0.067576 -12.5149 < 2.2e-16 ***
## sma1 -0.455487 0.064890 -7.0193 2.229e-12 ***
## sma2 -0.313974 0.062835 -4.9968 5.828e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Podemos observar que \(\theta_{1}\) no resulta significativa al \(5\%\). Por lo tanto ajustaremos nuevamente el modelo sin tenerla en cuenta y revisaremos si su ausencia genera algún cambio en la significancia de los demás parámetros.
fixed <- c(NA,NA,NA,NA,NA,0,NA,NA,NA)
fit2<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.038250 0.060448 -17.1758 < 2.2e-16 ***
## ar2 0.064321 0.141319 0.4551 0.6490
## ar3 0.225652 0.150466 1.4997 0.1337
## ar4 0.026229 0.124859 0.2101 0.8336
## ar5 -0.030537 0.075052 -0.4069 0.6841
## ma2 -0.765643 0.112103 -6.8298 8.503e-12 ***
## sma1 -0.448427 0.065652 -6.8303 8.472e-12 ***
## sma2 -0.302985 0.062386 -4.8566 1.194e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Después de ajustar el nuevo modelo podemos ver que \(\phi_{5}\) no resulta significativo mientras que \(\phi_{1}\) sigue siendo significativo al \(5\%\). Entonces ajustaremos un nuevo modelo omitiendo \(\phi_{5}\).
fixed <- c(NA,NA,NA,NA,0,0,NA,NA,NA)
fit2<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.039262 0.060385 -17.2105 < 2.2e-16 ***
## ar2 0.086090 0.132402 0.6502 0.5156
## ar3 0.253489 0.136654 1.8550 0.0636 .
## ar4 0.068419 0.072930 0.9381 0.3482
## ma2 -0.792228 0.090782 -8.7267 < 2.2e-16 ***
## sma1 -0.441893 0.063960 -6.9089 4.883e-12 ***
## sma2 -0.304659 0.062250 -4.8941 9.876e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,NA,NA,0,0,0,NA,NA,NA)
fit2<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.031370 0.059990 -17.1923 < 2.2e-16 ***
## ar2 0.038068 0.122224 0.3115 0.75545
## ar3 0.151554 0.079616 1.9036 0.05697 .
## ma2 -0.745319 0.081609 -9.1328 < 2.2e-16 ***
## sma1 -0.432562 0.062911 -6.8758 6.164e-12 ***
## sma2 -0.302525 0.062103 -4.8713 1.109e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,NA,0,0,0,0,NA,NA,NA)
fit2<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.045947 0.060512 -17.2849 < 2.2e-16 ***
## ar2 -0.157626 0.063384 -2.4868 0.01289 *
## ma2 -0.630294 0.080994 -7.7820 7.138e-15 ***
## sma1 -0.436232 0.062039 -7.0316 2.042e-12 ***
## sma2 -0.311758 0.061860 -5.0397 4.662e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Luego de refinar nuestro modelo podemos ver que solamente \(\phi_{1}, \phi_{2}\) y \(\theta_{2}\) son significativos al \(5\%\). Así, éste último será nuestro modelo seleccionado para la serie de tiempo Renta.
Estudio de Outliers
Nos centraremos ahora en saber si nuestra serie presenta outliers basándonos en los residuales, resultantes de nuestro modelo seleccionado anteriormente.
resi= residuals(fit2)
plot(resi)
coef= coefs2poly(fit2)
outliers= locate.outliers(resi,coef, cval = 4)
outliers## type ind coefhat tstat
## 1 AO 76 0.4496810 5.179170
## 2 AO 87 -0.3869329 -4.456434
## 3 AO 88 -0.8755647 -10.084167
## 4 AO 89 0.8270599 9.525522
## 5 AO 101 -0.4576588 -5.270939
## 6 AO 126 -0.4149417 -4.778729
## 7 AO 228 0.4983796 5.726962
## 8 AO 284 0.4770771 4.629161
## 10 LS 123 -0.2166965 -4.420095
## 16 TC 124 -0.3714931 -5.737153
## 17 TC 125 -0.2889362 -4.462179Podemos notar que existen 11 outliers en nuestra serie de tiempo, de los cuales la mayoría son outliers de cambio de nivel. Nuestro siguiente paso será obtener las intervenciones de éstos outliers y ajustar un nuevo modelo.
n=length(Renta)
xreg=outliers.effects(outliers,n)
xreg## AO76 AO87 AO88 AO89 AO101 AO126 AO228 AO284 LS123 TC124
## [1,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [2,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [3,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [4,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [5,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [6,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [7,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [8,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [9,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [10,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [11,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [12,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [13,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [14,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [15,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [16,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [17,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [18,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [19,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [20,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [21,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [22,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [23,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [24,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [25,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [26,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [27,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [28,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [29,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [30,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [31,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [32,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [33,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [34,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [35,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [36,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [37,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [38,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [39,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [40,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [41,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [42,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [43,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [44,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [45,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [46,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [47,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [48,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [49,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [50,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [51,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [52,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [53,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [54,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [55,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [56,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [57,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [58,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [59,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [60,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [61,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [62,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [63,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [64,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [65,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [66,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [67,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [68,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [69,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [70,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [71,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [72,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [73,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [74,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [75,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [76,] 1 0 0 0 0 0 0 0 0 0.000000e+00
## [77,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [78,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [79,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [80,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [81,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [82,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [83,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [84,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [85,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [86,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [87,] 0 1 0 0 0 0 0 0 0 0.000000e+00
## [88,] 0 0 1 0 0 0 0 0 0 0.000000e+00
## [89,] 0 0 0 1 0 0 0 0 0 0.000000e+00
## [90,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [91,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [92,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [93,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [94,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [95,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [96,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [97,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [98,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [99,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [100,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [101,] 0 0 0 0 1 0 0 0 0 0.000000e+00
## [102,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [103,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [104,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [105,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [106,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [107,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [108,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [109,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [110,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [111,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [112,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [113,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [114,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [115,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [116,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [117,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [118,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [119,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [120,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [121,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [122,] 0 0 0 0 0 0 0 0 0 0.000000e+00
## [123,] 0 0 0 0 0 0 0 0 1 0.000000e+00
## [124,] 0 0 0 0 0 0 0 0 1 1.000000e+00
## [125,] 0 0 0 0 0 0 0 0 1 7.000000e-01
## [126,] 0 0 0 0 0 1 0 0 1 4.900000e-01
## [127,] 0 0 0 0 0 0 0 0 1 3.430000e-01
## [128,] 0 0 0 0 0 0 0 0 1 2.401000e-01
## [129,] 0 0 0 0 0 0 0 0 1 1.680700e-01
## [130,] 0 0 0 0 0 0 0 0 1 1.176490e-01
## [131,] 0 0 0 0 0 0 0 0 1 8.235430e-02
## [132,] 0 0 0 0 0 0 0 0 1 5.764801e-02
## [133,] 0 0 0 0 0 0 0 0 1 4.035361e-02
## [134,] 0 0 0 0 0 0 0 0 1 2.824752e-02
## [135,] 0 0 0 0 0 0 0 0 1 1.977327e-02
## [136,] 0 0 0 0 0 0 0 0 1 1.384129e-02
## [137,] 0 0 0 0 0 0 0 0 1 9.688901e-03
## [138,] 0 0 0 0 0 0 0 0 1 6.782231e-03
## [139,] 0 0 0 0 0 0 0 0 1 4.747562e-03
## [140,] 0 0 0 0 0 0 0 0 1 3.323293e-03
## [141,] 0 0 0 0 0 0 0 0 1 2.326305e-03
## [142,] 0 0 0 0 0 0 0 0 1 1.628414e-03
## [143,] 0 0 0 0 0 0 0 0 1 1.139890e-03
## [144,] 0 0 0 0 0 0 0 0 1 7.979227e-04
## [145,] 0 0 0 0 0 0 0 0 1 5.585459e-04
## [146,] 0 0 0 0 0 0 0 0 1 3.909821e-04
## [147,] 0 0 0 0 0 0 0 0 1 2.736875e-04
## [148,] 0 0 0 0 0 0 0 0 1 1.915812e-04
## [149,] 0 0 0 0 0 0 0 0 1 1.341069e-04
## [150,] 0 0 0 0 0 0 0 0 1 9.387480e-05
## [151,] 0 0 0 0 0 0 0 0 1 6.571236e-05
## [152,] 0 0 0 0 0 0 0 0 1 4.599865e-05
## [153,] 0 0 0 0 0 0 0 0 1 3.219906e-05
## [154,] 0 0 0 0 0 0 0 0 1 2.253934e-05
## [155,] 0 0 0 0 0 0 0 0 1 1.577754e-05
## [156,] 0 0 0 0 0 0 0 0 1 1.104428e-05
## [157,] 0 0 0 0 0 0 0 0 1 7.730994e-06
## [158,] 0 0 0 0 0 0 0 0 1 5.411696e-06
## [159,] 0 0 0 0 0 0 0 0 1 3.788187e-06
## [160,] 0 0 0 0 0 0 0 0 1 2.651731e-06
## [161,] 0 0 0 0 0 0 0 0 1 1.856212e-06
## [162,] 0 0 0 0 0 0 0 0 1 1.299348e-06
## [163,] 0 0 0 0 0 0 0 0 1 9.095437e-07
## [164,] 0 0 0 0 0 0 0 0 1 6.366806e-07
## [165,] 0 0 0 0 0 0 0 0 1 4.456764e-07
## [166,] 0 0 0 0 0 0 0 0 1 3.119735e-07
## [167,] 0 0 0 0 0 0 0 0 1 2.183814e-07
## [168,] 0 0 0 0 0 0 0 0 1 1.528670e-07
## [169,] 0 0 0 0 0 0 0 0 1 1.070069e-07
## [170,] 0 0 0 0 0 0 0 0 1 7.490483e-08
## [171,] 0 0 0 0 0 0 0 0 1 5.243338e-08
## [172,] 0 0 0 0 0 0 0 0 1 3.670337e-08
## [173,] 0 0 0 0 0 0 0 0 1 2.569236e-08
## [174,] 0 0 0 0 0 0 0 0 1 1.798465e-08
## [175,] 0 0 0 0 0 0 0 0 1 1.258926e-08
## [176,] 0 0 0 0 0 0 0 0 1 8.812479e-09
## [177,] 0 0 0 0 0 0 0 0 1 6.168735e-09
## [178,] 0 0 0 0 0 0 0 0 1 4.318115e-09
## [179,] 0 0 0 0 0 0 0 0 1 3.022680e-09
## [180,] 0 0 0 0 0 0 0 0 1 2.115876e-09
## [181,] 0 0 0 0 0 0 0 0 1 1.481113e-09
## [182,] 0 0 0 0 0 0 0 0 1 1.036779e-09
## [183,] 0 0 0 0 0 0 0 0 1 7.257455e-10
## [184,] 0 0 0 0 0 0 0 0 1 5.080219e-10
## [185,] 0 0 0 0 0 0 0 0 1 3.556153e-10
## [186,] 0 0 0 0 0 0 0 0 1 2.489307e-10
## [187,] 0 0 0 0 0 0 0 0 1 1.742515e-10
## [188,] 0 0 0 0 0 0 0 0 1 1.219760e-10
## [189,] 0 0 0 0 0 0 0 0 1 8.538323e-11
## [190,] 0 0 0 0 0 0 0 0 1 5.976826e-11
## [191,] 0 0 0 0 0 0 0 0 1 4.183778e-11
## [192,] 0 0 0 0 0 0 0 0 1 2.928645e-11
## [193,] 0 0 0 0 0 0 0 0 1 2.050051e-11
## [194,] 0 0 0 0 0 0 0 0 1 1.435036e-11
## [195,] 0 0 0 0 0 0 0 0 1 1.004525e-11
## [196,] 0 0 0 0 0 0 0 0 1 7.031676e-12
## [197,] 0 0 0 0 0 0 0 0 1 4.922174e-12
## [198,] 0 0 0 0 0 0 0 0 1 3.445521e-12
## [199,] 0 0 0 0 0 0 0 0 1 2.411865e-12
## [200,] 0 0 0 0 0 0 0 0 1 1.688306e-12
## [201,] 0 0 0 0 0 0 0 0 1 1.181814e-12
## [202,] 0 0 0 0 0 0 0 0 1 8.272697e-13
## [203,] 0 0 0 0 0 0 0 0 1 5.790888e-13
## [204,] 0 0 0 0 0 0 0 0 1 4.053622e-13
## [205,] 0 0 0 0 0 0 0 0 1 2.837535e-13
## [206,] 0 0 0 0 0 0 0 0 1 1.986275e-13
## [207,] 0 0 0 0 0 0 0 0 1 1.390392e-13
## [208,] 0 0 0 0 0 0 0 0 1 9.732745e-14
## [209,] 0 0 0 0 0 0 0 0 1 6.812922e-14
## [210,] 0 0 0 0 0 0 0 0 1 4.769045e-14
## [211,] 0 0 0 0 0 0 0 0 1 3.338332e-14
## [212,] 0 0 0 0 0 0 0 0 1 2.336832e-14
## [213,] 0 0 0 0 0 0 0 0 1 1.635783e-14
## [214,] 0 0 0 0 0 0 0 0 1 1.145048e-14
## [215,] 0 0 0 0 0 0 0 0 1 8.015334e-15
## [216,] 0 0 0 0 0 0 0 0 1 5.610734e-15
## [217,] 0 0 0 0 0 0 0 0 1 3.927514e-15
## [218,] 0 0 0 0 0 0 0 0 1 2.749260e-15
## [219,] 0 0 0 0 0 0 0 0 1 1.924482e-15
## [220,] 0 0 0 0 0 0 0 0 1 1.347137e-15
## [221,] 0 0 0 0 0 0 0 0 1 9.429961e-16
## [222,] 0 0 0 0 0 0 0 0 1 6.600972e-16
## [223,] 0 0 0 0 0 0 0 0 1 4.620681e-16
## [224,] 0 0 0 0 0 0 0 0 1 3.234477e-16
## [225,] 0 0 0 0 0 0 0 0 1 2.264134e-16
## [226,] 0 0 0 0 0 0 0 0 1 1.584893e-16
## [227,] 0 0 0 0 0 0 0 0 1 1.109425e-16
## [228,] 0 0 0 0 0 0 1 0 1 7.765978e-17
## [229,] 0 0 0 0 0 0 0 0 1 5.436185e-17
## [230,] 0 0 0 0 0 0 0 0 1 3.805329e-17
## [231,] 0 0 0 0 0 0 0 0 1 2.663730e-17
## [232,] 0 0 0 0 0 0 0 0 1 1.864611e-17
## [233,] 0 0 0 0 0 0 0 0 1 1.305228e-17
## [234,] 0 0 0 0 0 0 0 0 1 9.136596e-18
## [235,] 0 0 0 0 0 0 0 0 1 6.395617e-18
## [236,] 0 0 0 0 0 0 0 0 1 4.476932e-18
## [237,] 0 0 0 0 0 0 0 0 1 3.133852e-18
## [238,] 0 0 0 0 0 0 0 0 1 2.193697e-18
## [239,] 0 0 0 0 0 0 0 0 1 1.535588e-18
## [240,] 0 0 0 0 0 0 0 0 1 1.074911e-18
## [241,] 0 0 0 0 0 0 0 0 1 7.524379e-19
## [242,] 0 0 0 0 0 0 0 0 1 5.267066e-19
## [243,] 0 0 0 0 0 0 0 0 1 3.686946e-19
## [244,] 0 0 0 0 0 0 0 0 1 2.580862e-19
## [245,] 0 0 0 0 0 0 0 0 1 1.806603e-19
## [246,] 0 0 0 0 0 0 0 0 1 1.264622e-19
## [247,] 0 0 0 0 0 0 0 0 1 8.852357e-20
## [248,] 0 0 0 0 0 0 0 0 1 6.196650e-20
## [249,] 0 0 0 0 0 0 0 0 1 4.337655e-20
## [250,] 0 0 0 0 0 0 0 0 1 3.036358e-20
## [251,] 0 0 0 0 0 0 0 0 1 2.125451e-20
## [252,] 0 0 0 0 0 0 0 0 1 1.487816e-20
## [253,] 0 0 0 0 0 0 0 0 1 1.041471e-20
## [254,] 0 0 0 0 0 0 0 0 1 7.290297e-21
## [255,] 0 0 0 0 0 0 0 0 1 5.103208e-21
## [256,] 0 0 0 0 0 0 0 0 1 3.572245e-21
## [257,] 0 0 0 0 0 0 0 0 1 2.500572e-21
## [258,] 0 0 0 0 0 0 0 0 1 1.750400e-21
## [259,] 0 0 0 0 0 0 0 0 1 1.225280e-21
## [260,] 0 0 0 0 0 0 0 0 1 8.576961e-22
## [261,] 0 0 0 0 0 0 0 0 1 6.003873e-22
## [262,] 0 0 0 0 0 0 0 0 1 4.202711e-22
## [263,] 0 0 0 0 0 0 0 0 1 2.941898e-22
## [264,] 0 0 0 0 0 0 0 0 1 2.059328e-22
## [265,] 0 0 0 0 0 0 0 0 1 1.441530e-22
## [266,] 0 0 0 0 0 0 0 0 1 1.009071e-22
## [267,] 0 0 0 0 0 0 0 0 1 7.063496e-23
## [268,] 0 0 0 0 0 0 0 0 1 4.944447e-23
## [269,] 0 0 0 0 0 0 0 0 1 3.461113e-23
## [270,] 0 0 0 0 0 0 0 0 1 2.422779e-23
## [271,] 0 0 0 0 0 0 0 0 1 1.695945e-23
## [272,] 0 0 0 0 0 0 0 0 1 1.187162e-23
## [273,] 0 0 0 0 0 0 0 0 1 8.310133e-24
## [274,] 0 0 0 0 0 0 0 0 1 5.817093e-24
## [275,] 0 0 0 0 0 0 0 0 1 4.071965e-24
## [276,] 0 0 0 0 0 0 0 0 1 2.850376e-24
## [277,] 0 0 0 0 0 0 0 0 1 1.995263e-24
## [278,] 0 0 0 0 0 0 0 0 1 1.396684e-24
## [279,] 0 0 0 0 0 0 0 0 1 9.776788e-25
## [280,] 0 0 0 0 0 0 0 0 1 6.843752e-25
## [281,] 0 0 0 0 0 0 0 0 1 4.790626e-25
## [282,] 0 0 0 0 0 0 0 0 1 3.353438e-25
## [283,] 0 0 0 0 0 0 0 0 1 2.347407e-25
## [284,] 0 0 0 0 0 0 0 1 1 1.643185e-25
## [285,] 0 0 0 0 0 0 0 0 1 1.150229e-25
## [286,] 0 0 0 0 0 0 0 0 1 8.051605e-26
## [287,] 0 0 0 0 0 0 0 0 1 5.636124e-26
## [288,] 0 0 0 0 0 0 0 0 1 3.945287e-26
## TC125
## [1,] 0.000000e+00
## [2,] 0.000000e+00
## [3,] 0.000000e+00
## [4,] 0.000000e+00
## [5,] 0.000000e+00
## [6,] 0.000000e+00
## [7,] 0.000000e+00
## [8,] 0.000000e+00
## [9,] 0.000000e+00
## [10,] 0.000000e+00
## [11,] 0.000000e+00
## [12,] 0.000000e+00
## [13,] 0.000000e+00
## [14,] 0.000000e+00
## [15,] 0.000000e+00
## [16,] 0.000000e+00
## [17,] 0.000000e+00
## [18,] 0.000000e+00
## [19,] 0.000000e+00
## [20,] 0.000000e+00
## [21,] 0.000000e+00
## [22,] 0.000000e+00
## [23,] 0.000000e+00
## [24,] 0.000000e+00
## [25,] 0.000000e+00
## [26,] 0.000000e+00
## [27,] 0.000000e+00
## [28,] 0.000000e+00
## [29,] 0.000000e+00
## [30,] 0.000000e+00
## [31,] 0.000000e+00
## [32,] 0.000000e+00
## [33,] 0.000000e+00
## [34,] 0.000000e+00
## [35,] 0.000000e+00
## [36,] 0.000000e+00
## [37,] 0.000000e+00
## [38,] 0.000000e+00
## [39,] 0.000000e+00
## [40,] 0.000000e+00
## [41,] 0.000000e+00
## [42,] 0.000000e+00
## [43,] 0.000000e+00
## [44,] 0.000000e+00
## [45,] 0.000000e+00
## [46,] 0.000000e+00
## [47,] 0.000000e+00
## [48,] 0.000000e+00
## [49,] 0.000000e+00
## [50,] 0.000000e+00
## [51,] 0.000000e+00
## [52,] 0.000000e+00
## [53,] 0.000000e+00
## [54,] 0.000000e+00
## [55,] 0.000000e+00
## [56,] 0.000000e+00
## [57,] 0.000000e+00
## [58,] 0.000000e+00
## [59,] 0.000000e+00
## [60,] 0.000000e+00
## [61,] 0.000000e+00
## [62,] 0.000000e+00
## [63,] 0.000000e+00
## [64,] 0.000000e+00
## [65,] 0.000000e+00
## [66,] 0.000000e+00
## [67,] 0.000000e+00
## [68,] 0.000000e+00
## [69,] 0.000000e+00
## [70,] 0.000000e+00
## [71,] 0.000000e+00
## [72,] 0.000000e+00
## [73,] 0.000000e+00
## [74,] 0.000000e+00
## [75,] 0.000000e+00
## [76,] 0.000000e+00
## [77,] 0.000000e+00
## [78,] 0.000000e+00
## [79,] 0.000000e+00
## [80,] 0.000000e+00
## [81,] 0.000000e+00
## [82,] 0.000000e+00
## [83,] 0.000000e+00
## [84,] 0.000000e+00
## [85,] 0.000000e+00
## [86,] 0.000000e+00
## [87,] 0.000000e+00
## [88,] 0.000000e+00
## [89,] 0.000000e+00
## [90,] 0.000000e+00
## [91,] 0.000000e+00
## [92,] 0.000000e+00
## [93,] 0.000000e+00
## [94,] 0.000000e+00
## [95,] 0.000000e+00
## [96,] 0.000000e+00
## [97,] 0.000000e+00
## [98,] 0.000000e+00
## [99,] 0.000000e+00
## [100,] 0.000000e+00
## [101,] 0.000000e+00
## [102,] 0.000000e+00
## [103,] 0.000000e+00
## [104,] 0.000000e+00
## [105,] 0.000000e+00
## [106,] 0.000000e+00
## [107,] 0.000000e+00
## [108,] 0.000000e+00
## [109,] 0.000000e+00
## [110,] 0.000000e+00
## [111,] 0.000000e+00
## [112,] 0.000000e+00
## [113,] 0.000000e+00
## [114,] 0.000000e+00
## [115,] 0.000000e+00
## [116,] 0.000000e+00
## [117,] 0.000000e+00
## [118,] 0.000000e+00
## [119,] 0.000000e+00
## [120,] 0.000000e+00
## [121,] 0.000000e+00
## [122,] 0.000000e+00
## [123,] 0.000000e+00
## [124,] 0.000000e+00
## [125,] 1.000000e+00
## [126,] 7.000000e-01
## [127,] 4.900000e-01
## [128,] 3.430000e-01
## [129,] 2.401000e-01
## [130,] 1.680700e-01
## [131,] 1.176490e-01
## [132,] 8.235430e-02
## [133,] 5.764801e-02
## [134,] 4.035361e-02
## [135,] 2.824752e-02
## [136,] 1.977327e-02
## [137,] 1.384129e-02
## [138,] 9.688901e-03
## [139,] 6.782231e-03
## [140,] 4.747562e-03
## [141,] 3.323293e-03
## [142,] 2.326305e-03
## [143,] 1.628414e-03
## [144,] 1.139890e-03
## [145,] 7.979227e-04
## [146,] 5.585459e-04
## [147,] 3.909821e-04
## [148,] 2.736875e-04
## [149,] 1.915812e-04
## [150,] 1.341069e-04
## [151,] 9.387480e-05
## [152,] 6.571236e-05
## [153,] 4.599865e-05
## [154,] 3.219906e-05
## [155,] 2.253934e-05
## [156,] 1.577754e-05
## [157,] 1.104428e-05
## [158,] 7.730994e-06
## [159,] 5.411696e-06
## [160,] 3.788187e-06
## [161,] 2.651731e-06
## [162,] 1.856212e-06
## [163,] 1.299348e-06
## [164,] 9.095437e-07
## [165,] 6.366806e-07
## [166,] 4.456764e-07
## [167,] 3.119735e-07
## [168,] 2.183814e-07
## [169,] 1.528670e-07
## [170,] 1.070069e-07
## [171,] 7.490483e-08
## [172,] 5.243338e-08
## [173,] 3.670337e-08
## [174,] 2.569236e-08
## [175,] 1.798465e-08
## [176,] 1.258926e-08
## [177,] 8.812479e-09
## [178,] 6.168735e-09
## [179,] 4.318115e-09
## [180,] 3.022680e-09
## [181,] 2.115876e-09
## [182,] 1.481113e-09
## [183,] 1.036779e-09
## [184,] 7.257455e-10
## [185,] 5.080219e-10
## [186,] 3.556153e-10
## [187,] 2.489307e-10
## [188,] 1.742515e-10
## [189,] 1.219760e-10
## [190,] 8.538323e-11
## [191,] 5.976826e-11
## [192,] 4.183778e-11
## [193,] 2.928645e-11
## [194,] 2.050051e-11
## [195,] 1.435036e-11
## [196,] 1.004525e-11
## [197,] 7.031676e-12
## [198,] 4.922174e-12
## [199,] 3.445521e-12
## [200,] 2.411865e-12
## [201,] 1.688306e-12
## [202,] 1.181814e-12
## [203,] 8.272697e-13
## [204,] 5.790888e-13
## [205,] 4.053622e-13
## [206,] 2.837535e-13
## [207,] 1.986275e-13
## [208,] 1.390392e-13
## [209,] 9.732745e-14
## [210,] 6.812922e-14
## [211,] 4.769045e-14
## [212,] 3.338332e-14
## [213,] 2.336832e-14
## [214,] 1.635783e-14
## [215,] 1.145048e-14
## [216,] 8.015334e-15
## [217,] 5.610734e-15
## [218,] 3.927514e-15
## [219,] 2.749260e-15
## [220,] 1.924482e-15
## [221,] 1.347137e-15
## [222,] 9.429961e-16
## [223,] 6.600972e-16
## [224,] 4.620681e-16
## [225,] 3.234477e-16
## [226,] 2.264134e-16
## [227,] 1.584893e-16
## [228,] 1.109425e-16
## [229,] 7.765978e-17
## [230,] 5.436185e-17
## [231,] 3.805329e-17
## [232,] 2.663730e-17
## [233,] 1.864611e-17
## [234,] 1.305228e-17
## [235,] 9.136596e-18
## [236,] 6.395617e-18
## [237,] 4.476932e-18
## [238,] 3.133852e-18
## [239,] 2.193697e-18
## [240,] 1.535588e-18
## [241,] 1.074911e-18
## [242,] 7.524379e-19
## [243,] 5.267066e-19
## [244,] 3.686946e-19
## [245,] 2.580862e-19
## [246,] 1.806603e-19
## [247,] 1.264622e-19
## [248,] 8.852357e-20
## [249,] 6.196650e-20
## [250,] 4.337655e-20
## [251,] 3.036358e-20
## [252,] 2.125451e-20
## [253,] 1.487816e-20
## [254,] 1.041471e-20
## [255,] 7.290297e-21
## [256,] 5.103208e-21
## [257,] 3.572245e-21
## [258,] 2.500572e-21
## [259,] 1.750400e-21
## [260,] 1.225280e-21
## [261,] 8.576961e-22
## [262,] 6.003873e-22
## [263,] 4.202711e-22
## [264,] 2.941898e-22
## [265,] 2.059328e-22
## [266,] 1.441530e-22
## [267,] 1.009071e-22
## [268,] 7.063496e-23
## [269,] 4.944447e-23
## [270,] 3.461113e-23
## [271,] 2.422779e-23
## [272,] 1.695945e-23
## [273,] 1.187162e-23
## [274,] 8.310133e-24
## [275,] 5.817093e-24
## [276,] 4.071965e-24
## [277,] 2.850376e-24
## [278,] 1.995263e-24
## [279,] 1.396684e-24
## [280,] 9.776788e-25
## [281,] 6.843752e-25
## [282,] 4.790626e-25
## [283,] 3.353438e-25
## [284,] 2.347407e-25
## [285,] 1.643185e-25
## [286,] 1.150229e-25
## [287,] 8.051605e-26
## [288,] 5.636124e-26Ajustando nuevamente el modelo teniendo en cuanta las intervenciones
fixed <- c(NA,NA,0,0,0,0,NA,NA,NA, rep(NA,11))
fit3<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.8456484 0.0634854 -13.3204 < 2.2e-16 ***
## ar2 -0.0096499 0.0771084 -0.1251 0.9004066
## ma2 -0.4343653 0.1218414 -3.5650 0.0003638 ***
## sma1 -0.0880958 0.0638684 -1.3793 0.1677916
## sma2 -0.3590600 0.0647687 -5.5437 2.961e-08 ***
## AO76 0.0475561 0.1008163 0.4717 0.6371331
## AO87 -0.1883147 0.0816254 -2.3071 0.0210515 *
## AO88 -0.9569520 0.1024599 -9.3398 < 2.2e-16 ***
## AO89 0.8763714 0.1076326 8.1422 3.880e-16 ***
## AO101 -0.0373872 0.1031476 -0.3625 0.7170059
## AO126 -0.4018806 0.0944562 -4.2547 2.094e-05 ***
## AO228 0.4779310 0.0778061 6.1426 8.119e-10 ***
## AO284 0.4150038 0.1111233 3.7346 0.0001880 ***
## LS123 -0.0869831 0.0746712 -1.1649 0.2440669
## TC124 -0.4252449 0.0947680 -4.4872 7.216e-06 ***
## TC125 0.2609030 0.0965382 2.7026 0.0068802 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ahora, procederemos a refinar el modelo.
fixed <- c(NA,NA,0,0,0,0,NA,NA,NA, rep(NA,8),0,NA,NA)
fit3<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.851329 0.063015 -13.5099 < 2.2e-16 ***
## ar2 -0.020632 0.076594 -0.2694 0.7876416
## ma2 -0.417569 0.120519 -3.4647 0.0005307 ***
## sma1 -0.091808 0.063523 -1.4453 0.1483813
## sma2 -0.356159 0.064326 -5.5368 3.081e-08 ***
## AO76 0.048057 0.100792 0.4768 0.6335076
## AO87 -0.178166 0.081264 -2.1924 0.0283480 *
## AO88 -0.956182 0.102349 -9.3424 < 2.2e-16 ***
## AO89 0.873832 0.107379 8.1379 4.023e-16 ***
## AO101 -0.035474 0.102984 -0.3445 0.7304956
## AO126 -0.409661 0.095461 -4.2914 1.776e-05 ***
## AO228 0.477541 0.078058 6.1177 9.491e-10 ***
## AO284 0.415021 0.111034 3.7378 0.0001857 ***
## TC124 -0.480430 0.083138 -5.7787 7.529e-09 ***
## TC125 0.254195 0.097464 2.6081 0.0091048 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,NA,0,0,0,0,NA,NA,NA, rep(NA,4),0,rep(NA,3),0,NA,NA)
fit3<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.849334 0.062733 -13.5389 < 2.2e-16 ***
## ar2 -0.018068 0.076431 -0.2364 0.8131275
## ma2 -0.420489 0.120352 -3.4938 0.0004761 ***
## sma1 -0.090862 0.063457 -1.4319 0.1521800
## sma2 -0.356382 0.064347 -5.5384 3.052e-08 ***
## AO76 0.047992 0.100876 0.4757 0.6342539
## AO87 -0.177578 0.081226 -2.1862 0.0287995 *
## AO88 -0.956161 0.102430 -9.3348 < 2.2e-16 ***
## AO89 0.896315 0.085323 10.5050 < 2.2e-16 ***
## AO126 -0.406221 0.094980 -4.2769 1.895e-05 ***
## AO228 0.477673 0.078060 6.1193 9.398e-10 ***
## AO284 0.414921 0.111121 3.7339 0.0001885 ***
## TC124 -0.480441 0.083140 -5.7787 7.528e-09 ***
## TC125 0.248602 0.096225 2.5835 0.0097791 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,NA,0,0,0,0,NA,NA,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA)
fit3<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.847981 0.062680 -13.5287 < 2.2e-16 ***
## ar2 -0.015185 0.075839 -0.2002 0.8413061
## ma2 -0.422400 0.119426 -3.5369 0.0004048 ***
## sma1 -0.091761 0.063383 -1.4477 0.1476914
## sma2 -0.355979 0.064323 -5.5342 3.126e-08 ***
## AO87 -0.177511 0.081299 -2.1834 0.0290039 *
## AO88 -0.986937 0.079445 -12.4228 < 2.2e-16 ***
## AO89 0.897025 0.085327 10.5128 < 2.2e-16 ***
## AO126 -0.406896 0.094954 -4.2852 1.826e-05 ***
## AO228 0.477917 0.078142 6.1160 9.593e-10 ***
## AO284 0.414604 0.111164 3.7297 0.0001917 ***
## TC124 -0.482344 0.083123 -5.8028 6.522e-09 ***
## TC125 0.250527 0.096138 2.6059 0.0091629 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,NA,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA)
fit3<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.831979 0.061616 -13.5027 < 2.2e-16 ***
## ar2 -0.010572 0.074758 -0.1414 0.8875376
## ma2 -0.408619 0.111614 -3.6610 0.0002512 ***
## sma2 -0.332948 0.063054 -5.2804 1.289e-07 ***
## AO87 -0.177877 0.077111 -2.3068 0.0210672 *
## AO88 -0.985076 0.075709 -13.0113 < 2.2e-16 ***
## AO89 0.895001 0.081234 11.0176 < 2.2e-16 ***
## AO126 -0.401058 0.090097 -4.4514 8.531e-06 ***
## AO228 0.474787 0.074035 6.4130 1.427e-10 ***
## AO284 0.406591 0.111339 3.6518 0.0002604 ***
## TC124 -0.485052 0.079194 -6.1248 9.078e-10 ***
## TC125 0.246825 0.091803 2.6886 0.0071742 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA)
fit3<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.826312 0.046556 -17.7488 < 2.2e-16 ***
## ma2 -0.417442 0.091732 -4.5507 5.347e-06 ***
## sma2 -0.332460 0.063005 -5.2768 1.315e-07 ***
## AO87 -0.177046 0.076934 -2.3013 0.0213769 *
## AO88 -0.984815 0.075761 -12.9990 < 2.2e-16 ***
## AO89 0.897958 0.078757 11.4016 < 2.2e-16 ***
## AO126 -0.402296 0.089422 -4.4989 6.832e-06 ***
## AO228 0.475489 0.073891 6.4350 1.235e-10 ***
## AO284 0.405768 0.111296 3.6458 0.0002665 ***
## TC124 -0.484730 0.079165 -6.1231 9.179e-10 ***
## TC125 0.247740 0.091287 2.7138 0.0066507 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ya que hemos terminado de refinar el modelo, es necesario hacer nuevamente la búsqueda de outliers en nuestra serie de tiempo. Para ello utilizaremos los residuales generados por nuestro nuevo modelo.
resi= residuals(fit3)
plot(resi)
coef= coefs2poly(fit3)
outliers2= locate.outliers(resi,coef, cval = 4)
outliers2## type ind coefhat tstat
## 1 AO 170 0.2758352 4.866408
## 2 AO 244 -0.2728730 -4.784422
## 3 LS 171 -0.1833856 -4.126106
## 5 LS 207 0.1931447 4.344098
## 6 LS 243 -0.1900203 -4.249576
## 8 TC 172 -0.2575829 -5.115417Nuevamente hemos encontrado 6 outliers en nuestra serie de tiempo por lo que nos encargaremos de implementar su intervención en nuestro modelo.
n=length(Renta)
xreg2=outliers.effects(outliers2,n)
xreg2## AO170 AO244 LS171 LS207 LS243 TC172
## [1,] 0 0 0 0 0 0.000000e+00
## [2,] 0 0 0 0 0 0.000000e+00
## [3,] 0 0 0 0 0 0.000000e+00
## [4,] 0 0 0 0 0 0.000000e+00
## [5,] 0 0 0 0 0 0.000000e+00
## [6,] 0 0 0 0 0 0.000000e+00
## [7,] 0 0 0 0 0 0.000000e+00
## [8,] 0 0 0 0 0 0.000000e+00
## [9,] 0 0 0 0 0 0.000000e+00
## [10,] 0 0 0 0 0 0.000000e+00
## [11,] 0 0 0 0 0 0.000000e+00
## [12,] 0 0 0 0 0 0.000000e+00
## [13,] 0 0 0 0 0 0.000000e+00
## [14,] 0 0 0 0 0 0.000000e+00
## [15,] 0 0 0 0 0 0.000000e+00
## [16,] 0 0 0 0 0 0.000000e+00
## [17,] 0 0 0 0 0 0.000000e+00
## [18,] 0 0 0 0 0 0.000000e+00
## [19,] 0 0 0 0 0 0.000000e+00
## [20,] 0 0 0 0 0 0.000000e+00
## [21,] 0 0 0 0 0 0.000000e+00
## [22,] 0 0 0 0 0 0.000000e+00
## [23,] 0 0 0 0 0 0.000000e+00
## [24,] 0 0 0 0 0 0.000000e+00
## [25,] 0 0 0 0 0 0.000000e+00
## [26,] 0 0 0 0 0 0.000000e+00
## [27,] 0 0 0 0 0 0.000000e+00
## [28,] 0 0 0 0 0 0.000000e+00
## [29,] 0 0 0 0 0 0.000000e+00
## [30,] 0 0 0 0 0 0.000000e+00
## [31,] 0 0 0 0 0 0.000000e+00
## [32,] 0 0 0 0 0 0.000000e+00
## [33,] 0 0 0 0 0 0.000000e+00
## [34,] 0 0 0 0 0 0.000000e+00
## [35,] 0 0 0 0 0 0.000000e+00
## [36,] 0 0 0 0 0 0.000000e+00
## [37,] 0 0 0 0 0 0.000000e+00
## [38,] 0 0 0 0 0 0.000000e+00
## [39,] 0 0 0 0 0 0.000000e+00
## [40,] 0 0 0 0 0 0.000000e+00
## [41,] 0 0 0 0 0 0.000000e+00
## [42,] 0 0 0 0 0 0.000000e+00
## [43,] 0 0 0 0 0 0.000000e+00
## [44,] 0 0 0 0 0 0.000000e+00
## [45,] 0 0 0 0 0 0.000000e+00
## [46,] 0 0 0 0 0 0.000000e+00
## [47,] 0 0 0 0 0 0.000000e+00
## [48,] 0 0 0 0 0 0.000000e+00
## [49,] 0 0 0 0 0 0.000000e+00
## [50,] 0 0 0 0 0 0.000000e+00
## [51,] 0 0 0 0 0 0.000000e+00
## [52,] 0 0 0 0 0 0.000000e+00
## [53,] 0 0 0 0 0 0.000000e+00
## [54,] 0 0 0 0 0 0.000000e+00
## [55,] 0 0 0 0 0 0.000000e+00
## [56,] 0 0 0 0 0 0.000000e+00
## [57,] 0 0 0 0 0 0.000000e+00
## [58,] 0 0 0 0 0 0.000000e+00
## [59,] 0 0 0 0 0 0.000000e+00
## [60,] 0 0 0 0 0 0.000000e+00
## [61,] 0 0 0 0 0 0.000000e+00
## [62,] 0 0 0 0 0 0.000000e+00
## [63,] 0 0 0 0 0 0.000000e+00
## [64,] 0 0 0 0 0 0.000000e+00
## [65,] 0 0 0 0 0 0.000000e+00
## [66,] 0 0 0 0 0 0.000000e+00
## [67,] 0 0 0 0 0 0.000000e+00
## [68,] 0 0 0 0 0 0.000000e+00
## [69,] 0 0 0 0 0 0.000000e+00
## [70,] 0 0 0 0 0 0.000000e+00
## [71,] 0 0 0 0 0 0.000000e+00
## [72,] 0 0 0 0 0 0.000000e+00
## [73,] 0 0 0 0 0 0.000000e+00
## [74,] 0 0 0 0 0 0.000000e+00
## [75,] 0 0 0 0 0 0.000000e+00
## [76,] 0 0 0 0 0 0.000000e+00
## [77,] 0 0 0 0 0 0.000000e+00
## [78,] 0 0 0 0 0 0.000000e+00
## [79,] 0 0 0 0 0 0.000000e+00
## [80,] 0 0 0 0 0 0.000000e+00
## [81,] 0 0 0 0 0 0.000000e+00
## [82,] 0 0 0 0 0 0.000000e+00
## [83,] 0 0 0 0 0 0.000000e+00
## [84,] 0 0 0 0 0 0.000000e+00
## [85,] 0 0 0 0 0 0.000000e+00
## [86,] 0 0 0 0 0 0.000000e+00
## [87,] 0 0 0 0 0 0.000000e+00
## [88,] 0 0 0 0 0 0.000000e+00
## [89,] 0 0 0 0 0 0.000000e+00
## [90,] 0 0 0 0 0 0.000000e+00
## [91,] 0 0 0 0 0 0.000000e+00
## [92,] 0 0 0 0 0 0.000000e+00
## [93,] 0 0 0 0 0 0.000000e+00
## [94,] 0 0 0 0 0 0.000000e+00
## [95,] 0 0 0 0 0 0.000000e+00
## [96,] 0 0 0 0 0 0.000000e+00
## [97,] 0 0 0 0 0 0.000000e+00
## [98,] 0 0 0 0 0 0.000000e+00
## [99,] 0 0 0 0 0 0.000000e+00
## [100,] 0 0 0 0 0 0.000000e+00
## [101,] 0 0 0 0 0 0.000000e+00
## [102,] 0 0 0 0 0 0.000000e+00
## [103,] 0 0 0 0 0 0.000000e+00
## [104,] 0 0 0 0 0 0.000000e+00
## [105,] 0 0 0 0 0 0.000000e+00
## [106,] 0 0 0 0 0 0.000000e+00
## [107,] 0 0 0 0 0 0.000000e+00
## [108,] 0 0 0 0 0 0.000000e+00
## [109,] 0 0 0 0 0 0.000000e+00
## [110,] 0 0 0 0 0 0.000000e+00
## [111,] 0 0 0 0 0 0.000000e+00
## [112,] 0 0 0 0 0 0.000000e+00
## [113,] 0 0 0 0 0 0.000000e+00
## [114,] 0 0 0 0 0 0.000000e+00
## [115,] 0 0 0 0 0 0.000000e+00
## [116,] 0 0 0 0 0 0.000000e+00
## [117,] 0 0 0 0 0 0.000000e+00
## [118,] 0 0 0 0 0 0.000000e+00
## [119,] 0 0 0 0 0 0.000000e+00
## [120,] 0 0 0 0 0 0.000000e+00
## [121,] 0 0 0 0 0 0.000000e+00
## [122,] 0 0 0 0 0 0.000000e+00
## [123,] 0 0 0 0 0 0.000000e+00
## [124,] 0 0 0 0 0 0.000000e+00
## [125,] 0 0 0 0 0 0.000000e+00
## [126,] 0 0 0 0 0 0.000000e+00
## [127,] 0 0 0 0 0 0.000000e+00
## [128,] 0 0 0 0 0 0.000000e+00
## [129,] 0 0 0 0 0 0.000000e+00
## [130,] 0 0 0 0 0 0.000000e+00
## [131,] 0 0 0 0 0 0.000000e+00
## [132,] 0 0 0 0 0 0.000000e+00
## [133,] 0 0 0 0 0 0.000000e+00
## [134,] 0 0 0 0 0 0.000000e+00
## [135,] 0 0 0 0 0 0.000000e+00
## [136,] 0 0 0 0 0 0.000000e+00
## [137,] 0 0 0 0 0 0.000000e+00
## [138,] 0 0 0 0 0 0.000000e+00
## [139,] 0 0 0 0 0 0.000000e+00
## [140,] 0 0 0 0 0 0.000000e+00
## [141,] 0 0 0 0 0 0.000000e+00
## [142,] 0 0 0 0 0 0.000000e+00
## [143,] 0 0 0 0 0 0.000000e+00
## [144,] 0 0 0 0 0 0.000000e+00
## [145,] 0 0 0 0 0 0.000000e+00
## [146,] 0 0 0 0 0 0.000000e+00
## [147,] 0 0 0 0 0 0.000000e+00
## [148,] 0 0 0 0 0 0.000000e+00
## [149,] 0 0 0 0 0 0.000000e+00
## [150,] 0 0 0 0 0 0.000000e+00
## [151,] 0 0 0 0 0 0.000000e+00
## [152,] 0 0 0 0 0 0.000000e+00
## [153,] 0 0 0 0 0 0.000000e+00
## [154,] 0 0 0 0 0 0.000000e+00
## [155,] 0 0 0 0 0 0.000000e+00
## [156,] 0 0 0 0 0 0.000000e+00
## [157,] 0 0 0 0 0 0.000000e+00
## [158,] 0 0 0 0 0 0.000000e+00
## [159,] 0 0 0 0 0 0.000000e+00
## [160,] 0 0 0 0 0 0.000000e+00
## [161,] 0 0 0 0 0 0.000000e+00
## [162,] 0 0 0 0 0 0.000000e+00
## [163,] 0 0 0 0 0 0.000000e+00
## [164,] 0 0 0 0 0 0.000000e+00
## [165,] 0 0 0 0 0 0.000000e+00
## [166,] 0 0 0 0 0 0.000000e+00
## [167,] 0 0 0 0 0 0.000000e+00
## [168,] 0 0 0 0 0 0.000000e+00
## [169,] 0 0 0 0 0 0.000000e+00
## [170,] 1 0 0 0 0 0.000000e+00
## [171,] 0 0 1 0 0 0.000000e+00
## [172,] 0 0 1 0 0 1.000000e+00
## [173,] 0 0 1 0 0 7.000000e-01
## [174,] 0 0 1 0 0 4.900000e-01
## [175,] 0 0 1 0 0 3.430000e-01
## [176,] 0 0 1 0 0 2.401000e-01
## [177,] 0 0 1 0 0 1.680700e-01
## [178,] 0 0 1 0 0 1.176490e-01
## [179,] 0 0 1 0 0 8.235430e-02
## [180,] 0 0 1 0 0 5.764801e-02
## [181,] 0 0 1 0 0 4.035361e-02
## [182,] 0 0 1 0 0 2.824752e-02
## [183,] 0 0 1 0 0 1.977327e-02
## [184,] 0 0 1 0 0 1.384129e-02
## [185,] 0 0 1 0 0 9.688901e-03
## [186,] 0 0 1 0 0 6.782231e-03
## [187,] 0 0 1 0 0 4.747562e-03
## [188,] 0 0 1 0 0 3.323293e-03
## [189,] 0 0 1 0 0 2.326305e-03
## [190,] 0 0 1 0 0 1.628414e-03
## [191,] 0 0 1 0 0 1.139890e-03
## [192,] 0 0 1 0 0 7.979227e-04
## [193,] 0 0 1 0 0 5.585459e-04
## [194,] 0 0 1 0 0 3.909821e-04
## [195,] 0 0 1 0 0 2.736875e-04
## [196,] 0 0 1 0 0 1.915812e-04
## [197,] 0 0 1 0 0 1.341069e-04
## [198,] 0 0 1 0 0 9.387480e-05
## [199,] 0 0 1 0 0 6.571236e-05
## [200,] 0 0 1 0 0 4.599865e-05
## [201,] 0 0 1 0 0 3.219906e-05
## [202,] 0 0 1 0 0 2.253934e-05
## [203,] 0 0 1 0 0 1.577754e-05
## [204,] 0 0 1 0 0 1.104428e-05
## [205,] 0 0 1 0 0 7.730994e-06
## [206,] 0 0 1 0 0 5.411696e-06
## [207,] 0 0 1 1 0 3.788187e-06
## [208,] 0 0 1 1 0 2.651731e-06
## [209,] 0 0 1 1 0 1.856212e-06
## [210,] 0 0 1 1 0 1.299348e-06
## [211,] 0 0 1 1 0 9.095437e-07
## [212,] 0 0 1 1 0 6.366806e-07
## [213,] 0 0 1 1 0 4.456764e-07
## [214,] 0 0 1 1 0 3.119735e-07
## [215,] 0 0 1 1 0 2.183814e-07
## [216,] 0 0 1 1 0 1.528670e-07
## [217,] 0 0 1 1 0 1.070069e-07
## [218,] 0 0 1 1 0 7.490483e-08
## [219,] 0 0 1 1 0 5.243338e-08
## [220,] 0 0 1 1 0 3.670337e-08
## [221,] 0 0 1 1 0 2.569236e-08
## [222,] 0 0 1 1 0 1.798465e-08
## [223,] 0 0 1 1 0 1.258926e-08
## [224,] 0 0 1 1 0 8.812479e-09
## [225,] 0 0 1 1 0 6.168735e-09
## [226,] 0 0 1 1 0 4.318115e-09
## [227,] 0 0 1 1 0 3.022680e-09
## [228,] 0 0 1 1 0 2.115876e-09
## [229,] 0 0 1 1 0 1.481113e-09
## [230,] 0 0 1 1 0 1.036779e-09
## [231,] 0 0 1 1 0 7.257455e-10
## [232,] 0 0 1 1 0 5.080219e-10
## [233,] 0 0 1 1 0 3.556153e-10
## [234,] 0 0 1 1 0 2.489307e-10
## [235,] 0 0 1 1 0 1.742515e-10
## [236,] 0 0 1 1 0 1.219760e-10
## [237,] 0 0 1 1 0 8.538323e-11
## [238,] 0 0 1 1 0 5.976826e-11
## [239,] 0 0 1 1 0 4.183778e-11
## [240,] 0 0 1 1 0 2.928645e-11
## [241,] 0 0 1 1 0 2.050051e-11
## [242,] 0 0 1 1 0 1.435036e-11
## [243,] 0 0 1 1 1 1.004525e-11
## [244,] 0 1 1 1 1 7.031676e-12
## [245,] 0 0 1 1 1 4.922174e-12
## [246,] 0 0 1 1 1 3.445521e-12
## [247,] 0 0 1 1 1 2.411865e-12
## [248,] 0 0 1 1 1 1.688306e-12
## [249,] 0 0 1 1 1 1.181814e-12
## [250,] 0 0 1 1 1 8.272697e-13
## [251,] 0 0 1 1 1 5.790888e-13
## [252,] 0 0 1 1 1 4.053622e-13
## [253,] 0 0 1 1 1 2.837535e-13
## [254,] 0 0 1 1 1 1.986275e-13
## [255,] 0 0 1 1 1 1.390392e-13
## [256,] 0 0 1 1 1 9.732745e-14
## [257,] 0 0 1 1 1 6.812922e-14
## [258,] 0 0 1 1 1 4.769045e-14
## [259,] 0 0 1 1 1 3.338332e-14
## [260,] 0 0 1 1 1 2.336832e-14
## [261,] 0 0 1 1 1 1.635783e-14
## [262,] 0 0 1 1 1 1.145048e-14
## [263,] 0 0 1 1 1 8.015334e-15
## [264,] 0 0 1 1 1 5.610734e-15
## [265,] 0 0 1 1 1 3.927514e-15
## [266,] 0 0 1 1 1 2.749260e-15
## [267,] 0 0 1 1 1 1.924482e-15
## [268,] 0 0 1 1 1 1.347137e-15
## [269,] 0 0 1 1 1 9.429961e-16
## [270,] 0 0 1 1 1 6.600972e-16
## [271,] 0 0 1 1 1 4.620681e-16
## [272,] 0 0 1 1 1 3.234477e-16
## [273,] 0 0 1 1 1 2.264134e-16
## [274,] 0 0 1 1 1 1.584893e-16
## [275,] 0 0 1 1 1 1.109425e-16
## [276,] 0 0 1 1 1 7.765978e-17
## [277,] 0 0 1 1 1 5.436185e-17
## [278,] 0 0 1 1 1 3.805329e-17
## [279,] 0 0 1 1 1 2.663730e-17
## [280,] 0 0 1 1 1 1.864611e-17
## [281,] 0 0 1 1 1 1.305228e-17
## [282,] 0 0 1 1 1 9.136596e-18
## [283,] 0 0 1 1 1 6.395617e-18
## [284,] 0 0 1 1 1 4.476932e-18
## [285,] 0 0 1 1 1 3.133852e-18
## [286,] 0 0 1 1 1 2.193697e-18
## [287,] 0 0 1 1 1 1.535588e-18
## [288,] 0 0 1 1 1 1.074911e-18for(i in 1:6){
xreg <- cbind(xreg,xreg2[,i])
}fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA,rep(NA,6))
fit4<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, fixed=fixed, xreg = xreg)
coeftest(fit4)## Warning in as.vector(est)/se: longitud de objeto mayor no es múltiplo de la
## longitud de uno menor## Warning in cbind(est, se, tval, pval): number of rows of result is not a
## multiple of vector length (arg 2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.895821 0.040285 -22.2371 < 2.2e-16 ***
## ar2 0.000000 0.101068 0.0000 1.000000
## ar3 0.000000 0.078677 0.0000 1.000000
## ar4 0.000000 0.067326 0.0000 1.000000
## ar5 0.000000 0.068368 0.0000 1.000000
## ma1 0.000000 0.070298 0.0000 1.000000
## ma2 -0.565730 0.082943 -6.8207 9.061e-12 ***
## sma1 0.000000 0.066235 0.0000 1.000000
## sma2 -0.438483 0.104043 -4.2144 2.504e-05 ***
## AO76 0.000000 0.073957 0.0000 1.000000
## AO87 -0.167889 0.084699 -1.9822 0.047458 *
## AO88 -0.986581 0.070723 -13.9499 < 2.2e-16 ***
## AO89 0.899295 0.070638 12.7311 < 2.2e-16 ***
## AO101 0.000000 0.060845 0.0000 1.000000
## AO126 -0.350544 0.049636 -7.0623 1.638e-12 ***
## AO228 0.514501 0.050678 10.1523 < 2.2e-16 ***
## AO284 0.422774 0.068643 6.1590 7.319e-10 ***
## LS123 0.000000 0.040285 0.0000 1.000000
## TC124 -0.419000 0.101068 -4.1457 3.388e-05 ***
## TC125 0.198238 0.078677 2.5197 0.011747 *
## 0.176760 0.067326 2.6255 0.008653 **
## -0.191823 0.068368 -2.8057 0.005020 **
## -0.011283 0.070298 -0.1605 0.872488
## 0.135228 0.082943 1.6304 0.103024
## -0.106088 0.066235 -1.6017 0.109221
## -0.182348 0.104043 -1.7526 0.079668 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA,rep(NA,5),0)
fit4<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, fixed=fixed, xreg = xreg)
coeftest(fit4)## Warning in as.vector(est)/se: longitud de objeto mayor no es múltiplo de la
## longitud de uno menor## Warning in cbind(est, se, tval, pval): number of rows of result is not a
## multiple of vector length (arg 2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.885675 0.042743 -20.7208 < 2.2e-16 ***
## ar2 0.000000 0.108156 0.0000 1.000000
## ar3 0.000000 0.075757 0.0000 1.000000
## ar4 0.000000 0.068253 0.0000 1.000000
## ar5 0.000000 0.068952 0.0000 1.000000
## ma1 0.000000 0.071273 0.0000 1.000000
## ma2 -0.532778 0.083016 -6.4178 1.383e-10 ***
## sma1 0.000000 0.066946 0.0000 1.000000
## sma2 -0.437629 0.104560 -4.1854 2.846e-05 ***
## AO76 0.000000 0.072239 0.0000 1.000000
## AO87 -0.175914 0.084478 -2.0824 0.037309 *
## AO88 -0.985953 0.074749 -13.1902 < 2.2e-16 ***
## AO89 0.895369 0.070929 12.6235 < 2.2e-16 ***
## AO101 0.000000 0.053435 0.0000 1.000000
## AO126 -0.353166 0.051911 -6.8032 1.023e-11 ***
## AO228 0.515722 0.052717 9.7828 < 2.2e-16 ***
## AO284 0.420154 0.042743 9.8297 < 2.2e-16 ***
## LS123 0.000000 0.108156 0.0000 1.000000
## TC124 -0.446897 0.075757 -5.8991 3.655e-09 ***
## TC125 0.205817 0.068253 3.0155 0.002565 **
## 0.180269 0.068952 2.6144 0.008938 **
## -0.203469 0.071273 -2.8548 0.004306 **
## -0.096301 0.083016 -1.1600 0.246033
## 0.137143 0.066946 2.0486 0.040505 *
## -0.102861 0.104560 -0.9838 0.325236
## 0.000000 0.072239 0.0000 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA,rep(NA,4),0,0)
fit4<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit4)## Warning in as.vector(est)/se: longitud de objeto mayor no es múltiplo de la
## longitud de uno menor## Warning in cbind(est, se, tval, pval): number of rows of result is not a
## multiple of vector length (arg 2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.868971 0.043933 -19.7794 < 2.2e-16 ***
## ar2 0.000000 0.105361 0.0000 1.000000
## ar3 0.000000 0.075534 0.0000 1.000000
## ar4 0.000000 0.068715 0.0000 1.000000
## ar5 0.000000 0.069103 0.0000 1.000000
## ma1 0.000000 0.071672 0.0000 1.000000
## ma2 -0.485692 0.083470 -5.8187 5.930e-09 ***
## sma1 0.000000 0.066366 0.0000 1.000000
## sma2 -0.429333 0.104335 -4.1150 3.873e-05 ***
## AO76 0.000000 0.072758 0.0000 1.000000
## AO87 -0.181108 0.085454 -2.1194 0.034060 *
## AO88 -0.984599 0.074941 -13.1383 < 2.2e-16 ***
## AO89 0.890003 0.067352 13.2141 < 2.2e-16 ***
## AO101 0.000000 0.056015 0.0000 1.000000
## AO126 -0.358520 0.052824 -6.7871 1.144e-11 ***
## AO228 0.500509 0.043933 11.3925 < 2.2e-16 ***
## AO284 0.423896 0.105361 4.0233 5.739e-05 ***
## LS123 0.000000 0.075534 0.0000 1.000000
## TC124 -0.451414 0.068715 -6.5694 5.053e-11 ***
## TC125 0.210836 0.069103 3.0511 0.002280 **
## 0.191066 0.071672 2.6658 0.007680 **
## -0.248904 0.083470 -2.9819 0.002864 **
## -0.089602 0.066366 -1.3501 0.176981
## 0.159050 0.104335 1.5244 0.127404
## 0.000000 0.072758 0.0000 1.000000
## 0.000000 0.085454 0.0000 1.000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA,rep(NA,3),0,0,0)
fit4<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit4)## Warning in as.vector(est)/se: longitud de objeto mayor no es múltiplo de la
## longitud de uno menor## Warning in cbind(est, se, tval, pval): number of rows of result is not a
## multiple of vector length (arg 2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.828234 0.046116 -17.9596 < 2.2e-16 ***
## ar2 0.000000 0.098073 0.0000 1.0000000
## ar3 0.000000 0.074717 0.0000 1.0000000
## ar4 0.000000 0.070377 0.0000 1.0000000
## ar5 0.000000 0.069376 0.0000 1.0000000
## ma1 0.000000 0.073011 0.0000 1.0000000
## ma2 -0.369603 0.084730 -4.3622 1.288e-05 ***
## sma1 0.000000 0.066883 0.0000 1.0000000
## sma2 -0.390112 0.102972 -3.7885 0.0001515 ***
## AO76 0.000000 0.074132 0.0000 1.0000000
## AO87 -0.191683 0.087092 -2.2009 0.0277409 *
## AO88 -0.980331 0.075936 -12.9099 < 2.2e-16 ***
## AO89 0.880448 0.067374 13.0680 < 2.2e-16 ***
## AO101 0.000000 0.062361 0.0000 1.0000000
## AO126 -0.379723 0.046116 -8.2340 < 2.2e-16 ***
## AO228 0.488333 0.098073 4.9793 6.382e-07 ***
## AO284 0.412778 0.074717 5.5246 3.303e-08 ***
## LS123 0.000000 0.070377 0.0000 1.0000000
## TC124 -0.466430 0.069376 -6.7232 1.778e-11 ***
## TC125 0.231622 0.073011 3.1724 0.0015116 **
## 0.218317 0.084730 2.5766 0.0099768 **
## -0.269072 0.066883 -4.0230 5.745e-05 ***
## -0.113878 0.102972 -1.1059 0.2687634
## 0.000000 0.074132 0.0000 1.0000000
## 0.000000 0.087092 0.0000 1.0000000
## 0.000000 0.075936 0.0000 1.0000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA,rep(NA,3),0,0,0)
fit4<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit4)## Warning in as.vector(est)/se: longitud de objeto mayor no es múltiplo de la
## longitud de uno menor## Warning in cbind(est, se, tval, pval): number of rows of result is not a
## multiple of vector length (arg 2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.828234 0.046116 -17.9596 < 2.2e-16 ***
## ar2 0.000000 0.098073 0.0000 1.0000000
## ar3 0.000000 0.074717 0.0000 1.0000000
## ar4 0.000000 0.070377 0.0000 1.0000000
## ar5 0.000000 0.069376 0.0000 1.0000000
## ma1 0.000000 0.073011 0.0000 1.0000000
## ma2 -0.369603 0.084730 -4.3622 1.288e-05 ***
## sma1 0.000000 0.066883 0.0000 1.0000000
## sma2 -0.390112 0.102972 -3.7885 0.0001515 ***
## AO76 0.000000 0.074132 0.0000 1.0000000
## AO87 -0.191683 0.087092 -2.2009 0.0277409 *
## AO88 -0.980331 0.075936 -12.9099 < 2.2e-16 ***
## AO89 0.880448 0.067374 13.0680 < 2.2e-16 ***
## AO101 0.000000 0.062361 0.0000 1.0000000
## AO126 -0.379723 0.046116 -8.2340 < 2.2e-16 ***
## AO228 0.488333 0.098073 4.9793 6.382e-07 ***
## AO284 0.412778 0.074717 5.5246 3.303e-08 ***
## LS123 0.000000 0.070377 0.0000 1.0000000
## TC124 -0.466430 0.069376 -6.7232 1.778e-11 ***
## TC125 0.231622 0.073011 3.1724 0.0015116 **
## 0.218317 0.084730 2.5766 0.0099768 **
## -0.269072 0.066883 -4.0230 5.745e-05 ***
## -0.113878 0.102972 -1.1059 0.2687634
## 0.000000 0.074132 0.0000 1.0000000
## 0.000000 0.087092 0.0000 1.0000000
## 0.000000 0.075936 0.0000 1.0000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixed <- c(NA,0,0,0,0,0,NA,0,NA,0, rep(NA,3),0,rep(NA,3),0,NA,NA,rep(NA,2),0,0,0,0)
fit4<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, include.drift = TRUE, fixed=fixed, xreg = xreg)## Warning in Arima(Renta, order = c(5, 1, 2), seasonal = list(order = c(0, : No
## drift term fitted as the order of difference is 2 or more.coeftest(fit4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.806020 0.046432 -17.3592 < 2.2e-16 ***
## ar2 0.000000 0.090880 0.0000 1.0000000
## ar3 0.000000 0.072012 0.0000 1.0000000
## ar4 0.000000 0.071063 0.0000 1.0000000
## ar5 0.000000 0.069224 0.0000 1.0000000
## ma1 0.000000 0.073451 0.0000 1.0000000
## ma2 -0.305712 0.084763 -3.6067 0.0003101 ***
## sma1 0.000000 0.067059 0.0000 1.0000000
## sma2 -0.365452 0.101808 -3.5896 0.0003312 ***
## AO76 0.000000 0.073967 0.0000 1.0000000
## AO87 -0.193986 0.087701 -2.2119 0.0269728 *
## AO88 -0.977227 0.067636 -14.4483 < 2.2e-16 ***
## AO89 0.876100 0.067342 13.0098 < 2.2e-16 ***
## AO101 0.000000 0.046432 0.0000 1.0000000
## AO126 -0.392797 0.090880 -4.3221 1.545e-05 ***
## AO228 0.487594 0.072012 6.7710 1.279e-11 ***
## AO284 0.405815 0.071063 5.7107 1.125e-08 ***
## LS123 0.000000 0.069224 0.0000 1.0000000
## TC124 -0.481958 0.073451 -6.5616 5.323e-11 ***
## TC125 0.241632 0.084763 2.8507 0.0043624 **
## 0.282760 0.067059 4.2166 2.480e-05 ***
## -0.271033 0.101808 -2.6622 0.0077633 **
## 0.000000 0.073967 0.0000 1.0000000
## 0.000000 0.087701 0.0000 1.0000000
## 0.000000 0.067636 0.0000 1.0000000
## 0.000000 0.067342 0.0000 1.0000000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Después de refinar nuestro modelo hacemos nuevamente la búsqueda de outliers en nuestra seris de tiempo.
resi= residuals(fit4)
plot(resi)
coef= coefs2poly(fit4)
outliers3= locate.outliers(resi,coef, cval = 4)
outliers3## [1] type ind coefhat tstat
## <0 rows> (or 0-length row.names)Finalmente, parece que hemos removido todos los outliers de nuestra serie. El siguiente paso será hacer el pronóstico de nuestro modelo basado en el último modelo construido fit4
Verificación de los supuestos
### Pronostico del modelo
#x11()
residuales <- fit4$residuals
plot(residuales)
acf(residuales,lag.max = 45)
pacf(residuales,lag.max = 45)
#Test de autocorrelaci?n
Box.test(residuales, lag = (length(residuales)/4), type = "Ljung-Box", fitdf = 2)##
## Box-Ljung test
##
## data: residuales
## X-squared = 64.95, df = 70, p-value = 0.6483######Análisis de Outliers
#Test de normalidad
jarque.bera.test(residuales)##
## Jarque Bera Test
##
## data: residuales
## X-squared = 35.429, df = 2, p-value = 2.027e-08###Estad?ticas CUSUM
res=residuales
cum=cumsum(res)/sd(res)
N=length(res)
cumq=cumsum(res^2)/sum(res^2)
Af=0.948 ###Cuantil del 95% para la estad?stica cusum
co=0.09821####Valor del cuantil aproximado para cusumsq para n/2
LS=Af*sqrt(N)+2*Af*c(1:length(res))/sqrt(N)
LI=-LS
LQS=co+(1:length(res))/N
LQI=-co+(1:length(res))/N
x11()
par(mfrow=c(2,1))
plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")
lines(LI,type="S",col="red")
#CUSUM Square
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")
lines(LQS,type="S",col="red")
lines(LQI,type="S",col="red")plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")
lines(LI,type="S",col="red")
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")
lines(LQS,type="S",col="red")
lines(LQI,type="S",col="red")
Rolling para el pronóstico un paso adelante
Por último, es importante conocer el Error de generalización del modelo SARIMA haciendo la predicción un paso adelante. Para ello, utilizamos el método Rolling dividiendo nuestro conjunto de datos en Entrenamiento(\(90\%\)) y Prueba(\(10\%\)).
h <- 1
train <- window(Renta,start=c(2000,01),end=c(2019,03))
test <- window(Renta,start=c(2019,04),end=c(2023,12))
n <- length(test) - h + 1
fixed <- c(NA,NA,0,0,0,0,NA,NA,NA)
fitmodelo<-Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, fixed=fixed)
fc <- ts(numeric(n), start=c(2019,04), end = c(2023,12), freq=12)
for(i in 1:n)
{
x <- window(Renta, end=c(2019.167) + (i-1)/12)
refit <- forecast::Arima(x, model=fitmodelo)
fc[i] <- forecast::forecast(refit, h=h)$mean[h]
}
dife=(test-fc)^2
ecm=(1/(length(test)))*sum(dife)
sqrt(ecm)## [1] 1679055df <- data.frame(
Time = seq(as.Date("2019-04-01"), as.Date("2023-12-01"), by = "month"),
Observed = test,
Predicted = fc
)library(ggplot2)
# Crear gráfico comparativo
ggplot(df, aes(x = Time)) +
# Línea para valores observados
geom_line(aes(y = Observed, color = "Observados"), size = 1.2) +
# Línea para valores predichos
geom_line(aes(y = Predicted, color = "Predichos"), size = 1.2) +
# Título y etiquetas
labs(
title = "Comparación entre Valores Observados y Predichos",
subtitle = "Modelo SARIMA",
x = "Tiempo",
y = "Valores"
) +
# Personalización de colores
scale_color_manual(
name = "",
values = c("Observados" = "#1f77b4", "Predichos" = "#ff7f0e")
) +
# Fondo más minimalista
theme_minimal(base_size = 15) +
# Mejorar el estilo del título
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "top"
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.
Comparación de modelos
Por último compararemos los cinco modelos que hemos ajustado a nuestra serie de tiempo Renta. Esto lo haremos teniendo en cuenta el Error de Generalización que obtuvimos para cada modelo haciendo la predicción un paso adelante. Para comparar nuestros modelos observemos la siguiente tabla:
| Modelo | Error de generalización |
|---|---|
| Suav.Exp | 861527.30 |
| Árboles | 2934669.31 |
| RN Multicapa | 2352524.54 |
| RNN | 8694792.69 |
| GRU | 10895734.99 |
| LSTM | 14666353.01 |
| SARIMA | 1679055 |
En base a estos resultados podemos concluir que el mejor modelo que presenta una mejor predicción un paso adelante para nuestra serie Renta es el modelo de Suavizamiento Exponencial. Así que, en base a este resultado, haremos el pronóstico \(h=12\) pasos adelante para nuestra series del impuesto sobre la renta.
Pronóstico con el mejor modelo
Rentats <- ts(Renta, frequency = 12, start = c(2000,01,01))
modelo_hw <- HoltWinters(Rentats, alpha = 0.1044525, beta = 0.04870827, gamma = 0.4508742)
# Realizando el pronóstico 12 pasos adelante
pronostico_hw <- forecast::forecast(modelo_hw, h = 12)
plot(pronostico_hw)
print(pronostico_hw)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2024 13609271 11987818 15230725 11129472 16089071
## Feb 2024 11954958 10323805 13586111 9460325 14449591
## Mar 2024 10559652 8917945 12201360 8048877 13070428
## Apr 2024 22106141 20452997 23759284 19577876 24634406
## May 2024 12239936 10574453 13905420 9692799 14787073
## Jun 2024 18856449 17177702 20535196 16289027 21423871
## Jul 2024 12961357 11268403 14654311 10372207 15550506
## Aug 2024 15267468 13559349 16975588 12655125 17879811
## Sep 2024 12403737 10679479 14127996 9766712 15040763
## Oct 2024 11829062 10087680 13570445 9165848 14492277
## Nov 2024 11395593 9636091 13155094 8704667 14086518
## Dec 2024 12182612 10403988 13961236 9462441 14902782fixed <- c(NA,NA,0,0,0,0,NA,NA,NA)
RentaSARIMA=forecast::Arima(Renta, order=c(5,1,2), seasonal = list(order=c(0,1,2), period=12), lambda = 0, fixed=fixed)
coeftest(RentaSARIMA)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -1.045947 0.060512 -17.2849 < 2.2e-16 ***
## ar2 -0.157626 0.063384 -2.4868 0.01289 *
## ma2 -0.630294 0.080994 -7.7820 7.138e-15 ***
## sma1 -0.436232 0.062039 -7.0316 2.042e-12 ***
## sma2 -0.311758 0.061860 -5.0397 4.662e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pronostico_Renta=forecast::forecast(RentaSARIMA,h=12,level=0.95)
plot(Pronostico_Renta)
Pronostico_Renta## Point Forecast Lo 95 Hi 95
## Jan 2024 13862236 9795372 19617589
## Feb 2024 11991512 8470372 16976391
## Mar 2024 8500256 5935350 12173562
## Apr 2024 27339804 19055054 39226593
## May 2024 11585144 8010773 16754384
## Jun 2024 24204138 16694501 35091812
## Jul 2024 12150293 8323836 17735766
## Aug 2024 16976975 11596697 24853430
## Sep 2024 12337159 8376758 18169977
## Oct 2024 11713618 7927992 17306885
## Nov 2024 10498089 7066613 15595855
## Dec 2024 11636830 7806543 17346451