ANALISIS DE INTERVENCION DE SERIES TEMPORALES UNIVARIANTES

Autor/a
Afiliación

Jaime Isaac

UNIVERSIDAD DE EL SALVADOR FACULTAD MULTIDISCIPLINARIA DE OCCIDENTE

ANALISIS DE INTERVENCION EN R SEGUNDO EJEMPLO.

Lectura de la base de datos.

library(haven)
Warning: package 'haven' was built under R version 4.3.3
store <- read_sav("store.sav")

pasando la data a serie de tiempo.

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
nortons <- ts(cbind(store[,1], store[,2]), 
              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_edmarts_1 <- nortons[,1]
nortons_edmarts_2 <- nortons[,2]

graficando las series de tiempo

library(TSstudio)
Warning: package 'TSstudio' was built under R version 4.3.3
ts_plot(nortons)
20002001200220032004200520062007394041424344454647
NortonsEdmartnortons
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)
2000200120022003200420052006200739404142434445
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
modelo_regresion <- lm(nortons_edmarts_1 ~ time(nortons_edmarts_1))

# Aplicar la prueba de Breusch-Pagan
resultado_bp <- bptest(modelo_regresion)

# 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

Interpretación:

  • 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
resultado <- adf.test(nortons_edmarts_1)
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)
2000200120022003200420052006200741424344454647
nortons_edmarts_2
tsdisplay(nortons_edmarts_1)

# Calcular el rango y la desviación estándar de la serie temporal
rango_serie <- max(nortons_edmarts_1) - min(nortons_edmarts_1)
desviacion_estandar <- sd(nortons_edmarts_1)

# Calcular el rango medio estandarizado
rango_medio_estandarizado <- rango_serie / desviacion_estandar
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
meses <- cycle(nortons_edmarts_1)

# Aplicar la prueba de Kruskal-Wallis
resultado_kruskal <- kruskal.test(nortons_edmarts_1~ factor(meses))

# 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

Interpretación del resultado:

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.

    Conclusión:

    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.

USANDO TSOLUTLIERS PARA DETECTAR OUTLIERS.

modelo=tsoutliers::tso(nortons_edmarts_1)
modelo1=tsoutliers::tso(nortons_edmarts_2)
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
modelonuevo=modelo$yadj
plot(modelo)

ts_plot(modelonuevo)
200020012002200320042005200620073838.53939.54040.54141.5
modelonuevo
library(zoo)
autoplot.zoo(modelonuevo,col='red')

library(writexl)
Warning: package 'writexl' was built under R version 4.3.3
modelonuevo=as.data.frame(modelonuevo)
write_xlsx(modelonuevo,'modelonuevo.xlsx')
library(forecast)
arima=auto.arima(nortons_edmarts_1,seasonal = TRUE)
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)
404244404244404244404244404244404244404244
nortons_edmarts_1 - Series (Y axis) vs. Lags (X axis)Lag 1Lag 2Lag 3Lag 4Lag 5Lag 6Lag 7Lag 8Lag 9Lag 10Lag 11Lag 12
library(TSstudio)
TSstudio::ts_cor(nortons_edmarts_1,lag=36)
00.510122436−0.200.20.40.6
Seasonal Lag 12Non-Seasonalnortons_edmarts_1 ACF and PACF PlotsLagACFPACF
# Crea la serie de tiempo con ceros hasta la observación 73 y 1 después
intervencion <- c(rep(0, 72), rep(1, 24))
intervencion=as.data.frame(intervencion)
# Combina las tres series de tiempo con cbind
nortons1 <- ts(cbind(store[,1], store[,2],intervencion), 
              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
nortons_edmarts_A <- nortons1[,1]
nortons_edmarts_B <- nortons1[,2]
intervencion=nortons1[,3]
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
modeloARI=Arima(nortons_edmarts_A,order = c(1,1,1),xreg = intervencion)
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
arch_test <- archTest(modeloARI$residuals)
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

Interpretación:

  • 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.

predicciones=forecast(modeloARI,1,xreg=intervencion)
summary(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)

predicciones$mean
          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
ultimo=nortons_edmarts_A[96]
ultimo
[1] 42.19159
pred=predicciones$mean
pred_rev=cumsum(pred)+ultimo
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
pred_rev_ts=ts(pred_rev,start = c(2008,1),frequency = 12)
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
store1 <- read_excel("store1.xlsx", sheet = "Hoja1", 
    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
m <- prophet(store1) 
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     
future <- make_future_dataframe(m, 
                    periods = 60) 
forecast <- predict(m, future) 
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
pred=forecast$yhat
pred=as.data.frame(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>.
Actual
Predicted
35
40
45
50
55
Jan 2000
Jan 2001
Jan 2002
Jan 2003
Jan 2004
Jan 2005
Jan 2006
Jan 2007
Jan 2008
#
#install.packages("CausalImpact")
# defining pre-intervention period
pre_period<-as.Date(c("2000-01-01", "2005-01-01"))
# defining post-intervention period
post_period<-as.Date(c("2006-02-01", "2007-01-01" ))

time_period = seq.Date(as.Date("2000-01-01"), as.Date("2007-01-01"), by = "1 month")

# Creating the time series elements
norton<- ts(store$Nortons)
edmart <- ts(store$Edmart)

# Merging the data
library(zoo)
serie<-zoo(cbind(norton, edmart), time_period)
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
impact <- CausalImpact(serie, pre_period, post_period)
# 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")