library(haven)
Warning: package 'haven' was built under R version 4.3.3
<- read_sav("store.sav") store
library(haven)
Warning: package 'haven' was built under R version 4.3.3
<- read_sav("store.sav") store
library(tsoutliers)
Warning: package 'tsoutliers' was built under R version 4.3.3
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
# Crear el objeto de la serie de tiempo
<- ts(cbind(store[,1], store[,2]),
nortons start = c(2000, 1), frequency = 12)
head(nortons)
Nortons Edmart
Jan 2000 40.00000 45.00000
Feb 2000 40.43144 45.66793
Mar 2000 38.90491 44.23079
Apr 2000 41.04949 45.74803
May 2000 40.79482 43.95047
Jun 2000 38.96505 45.09289
# Separar las columnas en dos series de tiempo
<- nortons[,1]
nortons_edmarts_1 <- nortons[,2] nortons_edmarts_2
library(TSstudio)
Warning: package 'TSstudio' was built under R version 4.3.3
ts_plot(nortons)
library(ggplot2)
Warning: package 'ggplot2' was built under R version 4.3.3
library(forecast)
Warning: package 'forecast' was built under R version 4.3.3
autoplot(nortons)
ts_plot(nortons_edmarts_1)
# Instalar y cargar el paquete lmtest
#install.packages("lmtest")
library(lmtest)
Warning: package 'lmtest' was built under R version 4.3.1
Loading required package: zoo
Warning: package 'zoo' was built under R version 4.3.1
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
# Ajustar un modelo de regresión simple para la serie temporal
<- lm(nortons_edmarts_1 ~ time(nortons_edmarts_1))
modelo_regresion
# Aplicar la prueba de Breusch-Pagan
<- bptest(modelo_regresion)
resultado_bp
# Mostrar los resultados de la prueba de Breusch-Pagan
print(resultado_bp)
studentized Breusch-Pagan test
data: modelo_regresion
BP = 9.84, df = 1, p-value = 0.001708
Hipótesis nula (H0): La varianza de los errores es constante (no hay heterocedasticidad).
Hipótesis alternativa (H1): La varianza de los errores no es constante (hay heterocedasticidad).
Si el valor p es menor a 0.05, puedes rechazar la hipótesis nula y concluir que hay heterocedasticidad.
library(tseries)
Warning: package 'tseries' was built under R version 4.3.3
<- adf.test(nortons_edmarts_1)
resultado print(resultado)
Augmented Dickey-Fuller Test
data: nortons_edmarts_1
Dickey-Fuller = -2.0823, Lag order = 4, p-value = 0.5424
alternative hypothesis: stationary
ts_plot(nortons_edmarts_2)
tsdisplay(nortons_edmarts_1)
# Calcular el rango y la desviación estándar de la serie temporal
<- max(nortons_edmarts_1) - min(nortons_edmarts_1)
rango_serie <- sd(nortons_edmarts_1)
desviacion_estandar
# Calcular el rango medio estandarizado
<- rango_serie / desviacion_estandar
rango_medio_estandarizado rango_medio_estandarizado
[1] 3.870219
Una serie temporal es estacionaria en varianza si su rango medio estandarizado se mantiene dentro de un rango razonable (generalmente cercano a 1). En general:
Si el rango medio estandarizado es cercano a 1, esto sugiere que la serie es estacionaria en varianza.
Si el rango medio estandarizado es significativamente mayor o menor que 1, esto sugiere que la serie no es estacionaria en varianza.
# Supongamos que tienes una serie temporal mensual
# Dividir la serie por meses o periodos
# Crear una variable de tiempo para agrupar por meses
<- cycle(nortons_edmarts_1)
meses
# Aplicar la prueba de Kruskal-Wallis
<- kruskal.test(nortons_edmarts_1~ factor(meses))
resultado_kruskal
# Imprimir el resultado
print(resultado_kruskal)
Kruskal-Wallis rank sum test
data: nortons_edmarts_1 by factor(meses)
Kruskal-Wallis chi-squared = 7.0032, df = 11, p-value = 0.7988
La hipótesis nula de la prueba de Kruskal-Wallis es que no hay diferencias significativas entre los grupos (en este caso, los diferentes meses), es decir, que la serie no presenta estacionalidad.
Valor p = 0.7988: Como el valor p es mucho mayor que el nivel de significancia típico (por ejemplo, 0.05), no tenemos evidencia suficiente para rechazar la hipótesis nula.
Dado que el valor p es elevado (0.7988), no podemos rechazar la hipótesis nula. Esto indica que no hay diferencias significativas entre los meses, y por lo tanto, el análisis no proporciona evidencia suficiente para concluir que la serie presenta estacionalidad.
=tsoutliers::tso(nortons_edmarts_1) modelo
=tsoutliers::tso(nortons_edmarts_2) modelo1
modelo
Series: nortons_edmarts_1
Regression with ARIMA(0,0,1) errors
Coefficients:
ma1 intercept LS73
-0.6637 39.9780 3.3459
s.e. 0.1044 0.0302 0.0635
sigma^2 = 0.5719: log likelihood = -108.17
AIC=224.33 AICc=224.77 BIC=234.59
Outliers:
type ind time coefhat tstat
1 LS 73 2006:01 3.346 52.71
modelo1
Series: nortons_edmarts_2
Regression with ARIMA(0,0,1) errors
Coefficients:
ma1 intercept LS74
-0.7401 44.9966 -2.4128
s.e. 0.0708 0.0226 0.0510
sigma^2 = 0.5311: log likelihood = -104.71
AIC=217.43 AICc=217.87 BIC=227.69
Outliers:
type ind time coefhat tstat
1 LS 74 2006:02 -2.413 -47.34
=modelo$yadj modelonuevo
plot(modelo)
ts_plot(modelonuevo)
library(zoo)
autoplot.zoo(modelonuevo,col='red')
library(writexl)
Warning: package 'writexl' was built under R version 4.3.3
=as.data.frame(modelonuevo)
modelonuevowrite_xlsx(modelonuevo,'modelonuevo.xlsx')
library(forecast)
=auto.arima(nortons_edmarts_1,seasonal = TRUE)
arima arima
Series: nortons_edmarts_1
ARIMA(1,1,1)
Coefficients:
ar1 ma1
-0.3589 -0.5630
s.e. 0.1250 0.1118
sigma^2 = 1.131: log likelihood = -140.09
AIC=286.18 AICc=286.45 BIC=293.85
ts_lags(nortons_edmarts_1)
library(TSstudio)
::ts_cor(nortons_edmarts_1,lag=36) TSstudio
# Crea la serie de tiempo con ceros hasta la observación 73 y 1 después
<- c(rep(0, 72), rep(1, 24))
intervencion =as.data.frame(intervencion)
intervencion# Combina las tres series de tiempo con cbind
<- ts(cbind(store[,1], store[,2],intervencion),
nortons1 start = c(2000, 1), frequency = 12)
# Imprime la nueva serie de tiempo
print(nortons1)
Nortons Edmart intervencion
Jan 2000 40.00000 45.00000 0
Feb 2000 40.43144 45.66793 0
Mar 2000 38.90491 44.23079 0
Apr 2000 41.04949 45.74803 0
May 2000 40.79482 43.95047 0
Jun 2000 38.96505 45.09289 0
Jul 2000 41.20665 45.64678 0
Aug 2000 39.21872 45.58908 0
Sep 2000 39.38229 44.52039 0
Oct 2000 39.98853 43.95227 0
Nov 2000 40.95767 45.09887 0
Dec 2000 38.52250 45.44805 0
Jan 2001 41.49853 45.57595 0
Feb 2001 39.41948 44.03433 0
Mar 2001 40.35938 44.80586 0
Apr 2001 39.08326 44.36263 0
May 2001 39.93370 45.55763 0
Jun 2001 39.69462 45.07795 0
Jul 2001 41.42150 45.28954 0
Aug 2001 39.38366 45.51096 0
Sep 2001 40.02709 43.49367 0
Oct 2001 40.68531 46.59243 0
Nov 2001 39.24654 44.79476 0
Dec 2001 41.18923 45.34552 0
Jan 2002 38.86948 43.48088 0
Feb 2002 39.30587 46.31593 0
Mar 2002 40.40389 45.77621 0
Apr 2002 40.82534 43.01576 0
May 2002 40.16791 45.76774 0
Jun 2002 38.90369 44.64137 0
Jul 2002 40.83819 44.81372 0
Aug 2002 39.09727 46.00716 0
Sep 2002 39.95560 45.42073 0
Oct 2002 41.00105 45.50654 0
Nov 2002 39.42877 43.96431 0
Dec 2002 39.12779 44.18657 0
Jan 2003 39.84701 46.30026 0
Feb 2003 41.59470 45.03862 0
Mar 2003 40.53879 44.19705 0
Apr 2003 39.59239 44.50918 0
May 2003 40.81107 45.77951 0
Jun 2003 40.06777 44.15344 0
Jul 2003 39.77795 46.38272 0
Aug 2003 39.16732 44.29628 0
Sep 2003 39.27702 45.95688 0
Oct 2003 40.84485 44.26073 0
Nov 2003 39.37070 44.58767 0
Dec 2003 41.16101 44.69048 0
Jan 2004 39.33956 45.62543 0
Feb 2004 39.15670 45.34470 0
Mar 2004 39.60650 43.57150 0
Apr 2004 40.07621 45.57024 0
May 2004 40.54080 46.35431 0
Jun 2004 38.85537 44.61505 0
Jul 2004 39.69195 44.55338 0
Aug 2004 40.28124 44.22279 0
Sep 2004 39.59780 46.72009 0
Oct 2004 40.74951 45.06759 0
Nov 2004 40.59315 44.13942 0
Dec 2004 39.70802 45.74584 0
Jan 2005 39.43060 44.26901 0
Feb 2005 40.18857 46.09534 0
Mar 2005 39.00475 43.96278 0
Apr 2005 41.63301 46.06950 0
May 2005 39.64983 45.00216 0
Jun 2005 39.25390 44.91395 0
Jul 2005 40.07919 44.60417 0
Aug 2005 41.28459 46.28756 0
Sep 2005 40.36234 43.68155 0
Oct 2005 39.28358 44.20174 0
Nov 2005 39.03326 45.58649 0
Dec 2005 40.33931 44.75545 0
Jan 2006 41.78165 43.02468 1
Feb 2006 42.89407 43.94695 1
Mar 2006 43.17365 42.91331 1
Apr 2006 44.45700 41.60731 1
May 2006 43.76628 42.58547 1
Jun 2006 43.06151 44.09846 1
Jul 2006 42.97146 40.52406 1
Aug 2006 44.40976 42.78494 1
Sep 2006 42.14661 42.04685 1
Oct 2006 44.70248 43.91654 1
Nov 2006 42.30960 42.32652 1
Dec 2006 44.20985 43.29694 1
Jan 2007 41.57983 42.05079 1
Feb 2007 44.14243 41.90897 1
Mar 2007 43.83169 43.39755 1
Apr 2007 42.18784 42.07177 1
May 2007 44.66558 43.52714 1
Jun 2007 43.47062 41.20087 1
Jul 2007 43.89885 42.48599 1
Aug 2007 41.25474 44.10092 1
Sep 2007 45.02576 42.92109 1
Oct 2007 42.58916 41.81701 1
Nov 2007 44.05089 42.49015 1
Dec 2007 42.19159 43.51575 1
# Separar las columnas en dos series de tiempo
<- nortons1[,1]
nortons_edmarts_A <- nortons1[,2]
nortons_edmarts_B =nortons1[,3] intervencion
auto.arima(nortons_edmarts_A)
Series: nortons_edmarts_A
ARIMA(1,1,1)
Coefficients:
ar1 ma1
-0.3589 -0.5630
s.e. 0.1250 0.1118
sigma^2 = 1.131: log likelihood = -140.09
AIC=286.18 AICc=286.45 BIC=293.85
=Arima(nortons_edmarts_A,order = c(1,1,1),xreg = intervencion) modeloARI
summary(modeloARI)
Series: nortons_edmarts_A
Regression with ARIMA(1,1,1) errors
Coefficients:
ar1 ma1 xreg
-0.4472 -1.0000 3.3319
s.e. 0.0918 0.0283 0.1286
sigma^2 = 0.6288: log likelihood = -113.99
AIC=235.99 AICc=236.43 BIC=246.2
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.03372176 0.7762539 0.6612063 -0.1198549 1.619161 0.4934638
ACF1
Training set -0.07337585
library(lmtest)
coeftest(modeloARI)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.447179 0.091788 -4.8719 1.105e-06 ***
ma1 -0.999995 0.028280 -35.3606 < 2.2e-16 ***
xreg 3.331943 0.128620 25.9054 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
checkresiduals(modeloARI)
Ljung-Box test
data: Residuals from Regression with ARIMA(1,1,1) errors
Q* = 16.767, df = 17, p-value = 0.4702
Model df: 2. Total lags used: 19
accuracy(modeloARI)
ME RMSE MAE MPE MAPE MASE
Training set -0.03372176 0.7762539 0.6612063 -0.1198549 1.619161 0.4934638
ACF1
Training set -0.07337585
library(tseries)
library(MTS)
Warning: package 'MTS' was built under R version 4.3.3
# Aplicar la prueba ARCH sobre los residuos
<- archTest(modeloARI$residuals) arch_test
Q(m) of squared series(LM test):
Test statistic: 4.40037 p-value: 0.9274837
Rank-based Test:
Test statistic: 6.375778 p-value: 0.7827662
# Imprimir el resultado
print(arch_test)
NULL
Hipótesis nula: La varianza de los residuos es constante (no hay heterocedasticidad, la serie es estacionaria en varianza).
Hipótesis alternativa: Hay heterocedasticidad, lo que indica que la varianza de los residuos no es constante.
Si el valor p es menor a un nivel de significancia (como 0.05), puedes rechazar la hipótesis nula y concluir que hay heterocedasticidad, es decir, la serie no es estacionaria en varianza.
Dado que el valor p es 0.9274837, mucho mayor que el nivel de significancia comúnmente usado (0.05), no podemos rechazar la hipótesis nula. Esto significa que no hay evidencia suficiente de heterocedasticidad en los residuos al cuadrado, lo que sugiere que la serie puede ser estacionaria en varianza.
=forecast(modeloARI,1,xreg=intervencion)
prediccionessummary(predicciones)
Forecast method: Regression with ARIMA(1,1,1) errors
Model Information:
Series: nortons_edmarts_A
Regression with ARIMA(1,1,1) errors
Coefficients:
ar1 ma1 xreg
-0.4472 -1.0000 3.3319
s.e. 0.0918 0.0283 0.1286
sigma^2 = 0.6288: log likelihood = -113.99
AIC=235.99 AICc=236.43 BIC=246.2
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.03372176 0.7762539 0.6612063 -0.1198549 1.619161 0.4934638
ACF1
Training set -0.07337585
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 2008 40.48321 39.46169 41.50473 38.92094 42.04548
Feb 2008 39.75719 38.64252 40.87185 38.05245 41.46192
Mar 2008 40.08185 38.94757 41.21613 38.34711 41.81659
Apr 2008 39.93667 38.79935 41.07398 38.19729 41.67604
May 2008 40.00159 38.86329 41.13989 38.26071 41.74247
Jun 2008 39.97256 38.83423 41.11089 38.23164 41.71348
Jul 2008 39.98554 38.84713 41.12395 38.24449 41.72659
Aug 2008 39.97973 38.84134 41.11813 38.23871 41.72075
Sep 2008 39.98233 38.84393 41.12073 38.24129 41.72337
Oct 2008 39.98117 38.84277 41.11957 38.24014 41.72220
Nov 2008 39.98169 38.84329 41.12009 38.24065 41.72272
Dec 2008 39.98146 38.84306 41.11986 38.24042 41.72249
Jan 2009 39.98156 38.84316 41.11996 38.24053 41.72259
Feb 2009 39.98151 38.84311 41.11991 38.24048 41.72255
Mar 2009 39.98153 38.84313 41.11994 38.24050 41.72257
Apr 2009 39.98153 38.84313 41.11993 38.24049 41.72256
May 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Jun 2009 39.98153 38.84313 41.11993 38.24049 41.72256
Jul 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Aug 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Sep 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Oct 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Nov 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Dec 2009 39.98153 38.84313 41.11993 38.24050 41.72256
Jan 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Feb 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Mar 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Apr 2010 39.98153 38.84313 41.11993 38.24050 41.72256
May 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Jun 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Jul 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Aug 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Sep 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Oct 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Nov 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Dec 2010 39.98153 38.84313 41.11993 38.24050 41.72256
Jan 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Feb 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Mar 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Apr 2011 39.98153 38.84313 41.11993 38.24050 41.72256
May 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Jun 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Jul 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Aug 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Sep 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Oct 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Nov 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Dec 2011 39.98153 38.84313 41.11993 38.24050 41.72256
Jan 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Feb 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Mar 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Apr 2012 39.98153 38.84313 41.11993 38.24050 41.72256
May 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Jun 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Jul 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Aug 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Sep 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Oct 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Nov 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Dec 2012 39.98153 38.84313 41.11993 38.24050 41.72256
Jan 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Feb 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Mar 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Apr 2013 39.98153 38.84313 41.11993 38.24050 41.72256
May 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Jun 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Jul 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Aug 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Sep 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Oct 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Nov 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Dec 2013 39.98153 38.84313 41.11993 38.24050 41.72256
Jan 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Feb 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Mar 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Apr 2014 43.31347 42.17507 44.45187 41.57244 45.05450
May 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Jun 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Jul 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Aug 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Sep 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Oct 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Nov 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Dec 2014 43.31347 42.17507 44.45187 41.57244 45.05450
Jan 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Feb 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Mar 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Apr 2015 43.31347 42.17507 44.45187 41.57244 45.05450
May 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Jun 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Jul 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Aug 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Sep 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Oct 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Nov 2015 43.31347 42.17507 44.45187 41.57244 45.05450
Dec 2015 43.31347 42.17507 44.45187 41.57244 45.05450
names(predicciones)
[1] "method" "model" "level" "mean" "lower" "upper"
[7] "x" "series" "fitted" "residuals"
library(ggplot2)
autoplot(predicciones)
$mean predicciones
Jan Feb Mar Apr May Jun Jul Aug
2008 40.48321 39.75719 40.08185 39.93667 40.00159 39.97256 39.98554 39.97973
2009 39.98156 39.98151 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153
2010 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153
2011 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153
2012 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153
2013 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153 39.98153
2014 43.31347 43.31347 43.31347 43.31347 43.31347 43.31347 43.31347 43.31347
2015 43.31347 43.31347 43.31347 43.31347 43.31347 43.31347 43.31347 43.31347
Sep Oct Nov Dec
2008 39.98233 39.98117 39.98169 39.98146
2009 39.98153 39.98153 39.98153 39.98153
2010 39.98153 39.98153 39.98153 39.98153
2011 39.98153 39.98153 39.98153 39.98153
2012 39.98153 39.98153 39.98153 39.98153
2013 39.98153 39.98153 39.98153 39.98153
2014 43.31347 43.31347 43.31347 43.31347
2015 43.31347 43.31347 43.31347 43.31347
=nortons_edmarts_A[96]
ultimo ultimo
[1] 42.19159
=predicciones$mean pred
=cumsum(pred)+ultimo pred_rev
pred_rev
[1] 82.6748 122.4320 162.5138 202.4505 242.4521 282.4247 322.4102
[8] 362.3899 402.3723 442.3534 482.3351 522.3166 562.2981 602.2796
[15] 642.2612 682.2427 722.2242 762.2058 802.1873 842.1688 882.1504
[22] 922.1319 962.1134 1002.0949 1042.0765 1082.0580 1122.0395 1162.0211
[29] 1202.0026 1241.9841 1281.9656 1321.9472 1361.9287 1401.9102 1441.8918
[36] 1481.8733 1521.8548 1561.8363 1601.8179 1641.7994 1681.7809 1721.7625
[43] 1761.7440 1801.7255 1841.7070 1881.6886 1921.6701 1961.6516 2001.6332
[50] 2041.6147 2081.5962 2121.5777 2161.5593 2201.5408 2241.5223 2281.5039
[57] 2321.4854 2361.4669 2401.4484 2441.4300 2481.4115 2521.3930 2561.3746
[64] 2601.3561 2641.3376 2681.3191 2721.3007 2761.2822 2801.2637 2841.2453
[71] 2881.2268 2921.2083 2964.5218 3007.8353 3051.1487 3094.4622 3137.7757
[78] 3181.0891 3224.4026 3267.7161 3311.0296 3354.3430 3397.6565 3440.9700
[85] 3484.2834 3527.5969 3570.9104 3614.2239 3657.5373 3700.8508 3744.1643
[92] 3787.4777 3830.7912 3874.1047 3917.4182 3960.7316
=ts(pred_rev,start = c(2008,1),frequency = 12) pred_rev_ts
ts.plot(nortons_edmarts_A,pred_rev_ts,lty=c(1,3),col=c(1,5))
library(prophet)
Warning: package 'prophet' was built under R version 4.3.3
Loading required package: Rcpp
Warning: package 'Rcpp' was built under R version 4.3.3
Loading required package: rlang
Warning: package 'rlang' was built under R version 4.3.3
library(readxl)
Warning: package 'readxl' was built under R version 4.3.3
<- read_excel("store1.xlsx", sheet = "Hoja1",
store1 col_types = c("date", "numeric"))
head(store1)
# A tibble: 6 × 2
ds y
<dttm> <dbl>
1 2000-01-01 00:00:00 40
2 2000-02-01 00:00:00 40.4
3 2000-03-01 00:00:00 38.9
4 2000-04-01 00:00:00 41.0
5 2000-05-01 00:00:00 40.8
6 2000-06-01 00:00:00 39.0
<- prophet(store1) m
Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
summary(m)
Length Class Mode
growth 1 -none- character
changepoints 25 POSIXct numeric
n.changepoints 1 -none- numeric
changepoint.range 1 -none- numeric
yearly.seasonality 1 -none- character
weekly.seasonality 1 -none- character
daily.seasonality 1 -none- character
holidays 0 -none- NULL
seasonality.mode 1 -none- character
seasonality.prior.scale 1 -none- numeric
changepoint.prior.scale 1 -none- numeric
holidays.prior.scale 1 -none- numeric
mcmc.samples 1 -none- numeric
interval.width 1 -none- numeric
uncertainty.samples 1 -none- numeric
specified.changepoints 1 -none- logical
start 1 POSIXct numeric
y.scale 1 -none- numeric
logistic.floor 1 -none- logical
t.scale 1 -none- numeric
changepoints.t 25 -none- numeric
seasonalities 1 -none- list
extra_regressors 0 -none- list
country_holidays 0 -none- NULL
stan.fit 4 -none- list
params 6 -none- list
history 5 tbl_df list
history.dates 96 POSIXct numeric
train.holiday.names 0 -none- NULL
train.component.cols 3 data.frame list
component.modes 2 -none- list
fit.kwargs 0 -none- list
<- make_future_dataframe(m,
future periods = 60)
<- predict(m, future)
forecast tail(forecast[c('ds', 'yhat',
'yhat_lower', 'yhat_upper')])
ds yhat yhat_lower yhat_upper
151 2008-01-25 47.58837 46.42289 48.78611
152 2008-01-26 46.97188 45.81183 48.17889
153 2008-01-27 46.40408 45.20464 47.56081
154 2008-01-28 45.89957 44.63341 46.91735
155 2008-01-29 45.47127 44.26718 46.73704
156 2008-01-30 45.13008 43.90960 46.31818
=forecast$yhat pred
=as.data.frame(pred) pred
pred
pred
1 39.28955
2 40.30732
3 39.56172
4 41.08064
5 40.43544
6 39.06763
7 40.62451
8 40.16971
9 39.38105
10 40.11665
11 39.84674
12 39.44341
13 40.13727
14 40.24887
15 39.89621
16 40.63276
17 40.49985
18 39.34892
19 40.47685
20 39.79931
21 39.64012
22 40.19795
23 39.66565
24 39.61994
25 39.83261
26 40.23885
27 40.22547
28 40.18162
29 40.57170
30 39.62817
31 40.32147
32 39.43394
33 39.90561
34 40.27476
35 39.48384
36 39.79439
37 39.52709
38 40.23549
39 40.54920
40 39.72760
41 40.65093
42 39.90507
43 40.15855
44 39.07402
45 40.17733
46 40.34704
47 39.30173
48 39.96714
49 39.22168
50 40.23977
51 39.49447
52 41.01371
53 40.36916
54 39.00201
55 40.55954
56 40.16806
57 39.44273
58 40.23961
59 40.03325
60 39.69141
61 40.44880
62 40.62409
63 40.32896
64 41.12921
65 41.11247
66 40.08158
67 41.32568
68 40.76823
69 40.72913
70 41.40318
71 40.99097
72 41.06148
73 41.39424
74 41.92058
75 42.01567
76 42.09191
77 42.59788
78 41.77411
79 42.58330
80 41.81553
81 42.40695
82 42.89200
83 42.22084
84 42.64728
85 42.49974
86 43.32790
87 43.74978
88 43.04793
89 44.08716
90 43.46106
91 43.83043
92 42.86562
93 44.08866
94 44.37423
95 43.44848
96 44.22960
97 43.46317
98 42.66683
99 41.84782
100 41.01423
101 40.17504
102 39.33996
103 38.51940
104 37.72424
105 36.96577
106 36.25542
107 35.60460
108 35.02447
109 34.52569
110 34.11824
111 33.81110
112 33.61211
113 33.52768
114 33.56268
115 33.72018
116 34.00138
117 34.40548
118 34.92959
119 35.56879
120 36.31605
121 37.16239
122 38.09691
123 39.10703
124 40.17863
125 41.29632
126 42.44371
127 43.60370
128 44.75881
129 45.89154
130 46.98469
131 48.02171
132 48.98707
133 49.86656
134 50.64760
135 51.31951
136 51.87375
137 52.30409
138 52.60679
139 52.78065
140 52.82707
141 52.75002
142 52.55596
143 52.25373
144 51.85430
145 51.37061
146 50.81724
147 50.21010
148 49.56609
149 48.90269
150 48.23762
151 47.58837
152 46.97188
153 46.40408
154 45.89957
155 45.47127
156 45.13008
plot(m,forecast)
prophet_plot_components(m, forecast)
dyplot.prophet(m, forecast)
Warning: `select_()` was deprecated in dplyr 0.7.0.
ℹ Please use `select()` instead.
ℹ The deprecated feature was likely used in the prophet package.
Please report the issue at <https://github.com/facebook/prophet/issues>.
#
#install.packages("CausalImpact")
# defining pre-intervention period
<-as.Date(c("2000-01-01", "2005-01-01"))
pre_period# defining post-intervention period
<-as.Date(c("2006-02-01", "2007-01-01" ))
post_period
= seq.Date(as.Date("2000-01-01"), as.Date("2007-01-01"), by = "1 month")
time_period
# Creating the time series elements
<- ts(store$Nortons)
norton<- ts(store$Edmart)
edmart
# Merging the data
library(zoo)
<-zoo(cbind(norton, edmart), time_period) serie
library(CausalImpact)
Warning: package 'CausalImpact' was built under R version 4.3.3
Loading required package: bsts
Warning: package 'bsts' was built under R version 4.3.3
Loading required package: BoomSpikeSlab
Warning: package 'BoomSpikeSlab' was built under R version 4.3.3
Loading required package: Boom
Warning: package 'Boom' was built under R version 4.3.3
Attaching package: 'Boom'
The following object is masked from 'package:stats':
rWishart
Attaching package: 'BoomSpikeSlab'
The following object is masked from 'package:stats':
knots
Loading required package: xts
Warning: package 'xts' was built under R version 4.3.3
Attaching package: 'bsts'
The following object is masked from 'package:BoomSpikeSlab':
SuggestBurn
<- CausalImpact(serie, pre_period, post_period)
impact # Summary
summary(impact)
Posterior inference {CausalImpact}
Average Cumulative
Actual 43 520
Prediction (s.d.) 40 (0.31) 477 (3.71)
95% CI [39, 40] [469, 484]
Absolute effect (s.d.) 3.6 (0.31) 43.0 (3.71)
95% CI [3, 4.2] [36, 50.8]
Relative effect (s.d.) 9.1% (0.85%) 9.1% (0.85%)
95% CI [7.5%, 11%] [7.5%, 11%]
Posterior tail-area probability p: 0.00102
Posterior prob. of a causal effect: 99.89827%
For more details, type: summary(impact, "report")
summary(impact,"report")
Analysis report {CausalImpact}
During the post-intervention period, the response variable had an average value of approx. 43.31. By contrast, in the absence of an intervention, we would have expected an average response of 39.72. The 95% interval of this counterfactual prediction is [39.07, 40.30]. Subtracting this prediction from the observed response yields an estimate of the causal effect the intervention had on the response variable. This effect is 3.59 with a 95% interval of [3.01, 4.24]. For a discussion of the significance of this effect, see below.
Summing up the individual data points during the post-intervention period (which can only sometimes be meaningfully interpreted), the response variable had an overall value of 519.68. By contrast, had the intervention not taken place, we would have expected a sum of 476.65. The 95% interval of this prediction is [468.86, 483.56].
The above results are given in terms of absolute numbers. In relative terms, the response variable showed an increase of +9%. The 95% interval of this percentage is [+7%, +11%].
This means that the positive effect observed during the intervention period is statistically significant and unlikely to be due to random fluctuations. It should be noted, however, that the question of whether this increase also bears substantive significance can only be answered by comparing the absolute effect (3.59) to the original goal of the underlying intervention.
The probability of obtaining this effect by chance is very small (Bayesian one-sided tail-area probability p = 0.001). This means the causal effect can be considered statistically significant.
plot(impact)+scale_x_date(date_labels = "%b-%y", date_breaks = "1 year")