MODELOS DE DESCOMPOSICION

• Es importante saber que no debemos estimar un componente en presencia del otro (dentencia y estacionalidad)
• Cuando en una data hay mes y año es porque es estacional
• t es un consecutivo: 1,2,3,4,5,6…
• Zt= el valor en ese t
• Con la funcion decompose el R calcula el componente estacional como medias moviles, no usa regresion

Tendencia

Metodo de regresion

Cargando liberias

library("readxl")

Cargando datos

serie1=read_excel('serie1.xlsx')
str(serie1)

## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
##  $ t : num [1:60] 1 2 3 4 5 6 7 8 9 10 ...
##  $ t2: num [1:60] 1 4 9 16 25 36 49 64 81 100 ...
##  $ Zt: num [1:60] 390 378 358 364 392 ...

Graficando la serie

t=serie1$t
t2=serie1$t2
Zt=serie1$Zt
plot(t,Zt)

#Se observa una tendencia cuadratica

Sólo con fines ilustrativos y comparativos ejecutaremos una regresión lineal y otra cuadrática

El modelo lineal:

regre1=lm(Zt~t,data=serie1)
summary(regre1) 

## 
## Call:
## lm(formula = Zt ~ t, data = serie1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -118.088  -49.713   -1.655   43.711  115.379 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 287.35571   16.33960  17.586   <2e-16 ***
## t            -0.03336    0.46586  -0.072    0.943    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 62.49 on 58 degrees of freedom
## Multiple R-squared:  8.838e-05,  Adjusted R-squared:  -0.01715 
## F-statistic: 0.005127 on 1 and 58 DF,  p-value: 0.9432

#grafica
plot(t,Zt)
abline(reg=regre1, col="blue")

El modelo cuadratico (regresion):

regre2=lm(Zt~t+t2,data=serie1)
summary(regre2)

## 
## Call:
## lm(formula = Zt ~ t + t2, data = serie1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -55.428 -18.391   1.646  20.208  42.836 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 420.04498    9.97564   42.11   <2e-16 ***
## t           -12.87425    0.75458  -17.06   <2e-16 ***
## t2            0.21051    0.01199   17.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.9 on 57 degrees of freedom
## Multiple R-squared:  0.844,  Adjusted R-squared:  0.8385 
## F-statistic: 154.1 on 2 and 57 DF,  p-value: < 2.2e-16

#grafica
min_xlim=min(t)
max_xlim=max(t)+6 #prediccion de 6 puntos al futuro

datan=data.frame(t=seq(min_xlim,max_xlim),
                 t2=(seq(min_xlim,max_xlim)**2))
datan$y=predict(regre2,newdata = datan)
max_ylim=max(datan$y)
plot(t,Zt,xlim = c(min_xlim,max_xlim),ylim = c(min_xlim,max_ylim))
lines(x=datan$t,y=datan$y,col="red")

Metodo de suavización

• El filtro es una funcion o calculo: media, media ponderada, etc que requiere un n
• Para un filtro que use media poderada,la condicion es que el peso total de las observaciones sumen 1 pero no se divide entre la cantidad de elementos
• El numero de observaciones que intervienen en el filtro (calculo) es 2n-1
• Se utiliza lo que llamamos un filtro centrado. n hacia atras y n hacia adelante de toda la serie
• ejemplo tengo un vector serie=c(389, 377, 358, 363, 391, 330, 329,325..) y hago un filtro n=2, entoces tomare los 5 primeros datos (2*2)+1 = 2 y el resultado sera para estimar el dato del centro, que para el ejemplo sera: 358 y de esa manera sucesivamente
• Se pierden n hacia adelante de y n hacia atras de toda la serie
• Al utilizar promedio movil ser el promedio de los valores
• La varianza de X promedio es igual al X/n. Por lo tanto los calculo de filtros reducen la varibilidad de la serie original
• Esta opcion no permite pronosticar

library(stats)
n = 5
filtro1=filter(Zt, rep(1 / n, n), sides = 2)    # filtro1 (media moviles o pronedio movil simple). sides = 2 es centrado
filtro1

## Time Series:
## Start = 1 
## End = 60 
## Frequency = 1 
##  [1]     NA     NA 376.18 364.40 354.72 353.78 345.70 337.32 336.70 329.14
## [11] 318.62 315.18 301.34 292.70 291.40 276.56 266.84 265.70 262.56 256.70
## [21] 257.82 247.70 240.18 226.74 212.48 201.00 193.28 186.18 185.44 187.56
## [31] 204.68 212.26 218.98 215.58 223.18 224.08 226.08 235.90 252.40 260.94
## [41] 263.70 280.28 276.24 280.52 289.92 293.20 288.70 297.20 298.92 297.06
## [51] 303.94 318.20 329.16 338.90 343.42 357.68 364.80 371.68     NA     NA

library(stats)
#Filtro de 7 términos centrado con ponderaciones
filtro2 = filter(Zt, c(0.05,0.20,0.5,0.20,0.05), sides = 2)  # filtro2 (media ponderada o pronedio movil ponderado) sides = 2 es centrado
filtro2

## Time Series:
## Start = 1 
## End = 60 
## Frequency = 1 
##  [1]      NA      NA 366.415 367.220 369.075 345.455 337.285 341.265 335.125
## [10] 337.300 323.315 304.290 300.350 298.115 287.665 282.005 271.670 249.785
## [19] 256.725 266.780 265.005 252.185 233.940 220.200 221.135 204.225 187.805
## [28] 178.125 178.570 192.585 199.745 209.605 232.120 221.040 219.745 207.685
## [37] 228.285 245.575 242.785 265.155 268.485 275.495 280.710 286.160 278.955
## [46] 289.510 304.420 297.005 289.470 297.780 301.045 314.060 329.775 344.930
## [55] 352.030 348.005 353.880 381.220      NA      NA

Graficando

series=data.frame(Zt,filtro1,filtro2)
#stSeries=ts(series,freq=12,start=c(1981,1))
stSeries=ts(series,freq=1,start=1981)
plot(stSeries)

## Metodo de diferencias

Este metodo no es para obtner una funcion de tendencia, sino para remover tendencia. Si tengo un tendencia polinomica de grado 2 y se hace diferencia de t - t-1, t-1 -t2…en algun momento llegara a una constante. Siempre se perdera el primer dato. Baja como escalera

• Lo cubico en una primera diferencia se vuelve cuadratico, en una segunda diferencia se vuelve lineal y el una tercera diferencia se vuelve constante (revisar graficamente). La diferenciacion va haciendo que el grado de polinomio disminuya

• El resumen con este metodo podemos remover tendencia

• Las diferencias introducen correlaciones entre los residuos ya no se volveria un ruido blanco

d1=diff(Zt) #Una primera diferncia 
d1

##  [1] -11.9 -19.7   5.6  27.9 -60.8  -1.4  24.0 -30.1  26.4 -22.0 -36.1   9.2
## [13]   5.3 -25.6   4.0   0.6 -58.5  30.9  17.3  -6.0 -13.0 -23.6 -25.3  30.3
## [25] -35.6 -17.1  -9.7  -6.5  33.4  -3.8  -2.8  65.3 -54.2  29.1 -54.4  52.2
## [37]  31.8 -48.7  68.2 -21.0  12.4   2.9  20.4 -34.9  20.6  38.0 -27.7 -18.5
## [49]  30.1 -13.3  20.1  16.0  18.4  13.6 -19.4  -6.0  64.7 -17.3  12.4

d2=diff(Zt, differences = 2) #Una segunda diferencia 
d2

##  [1]   -7.8   25.3   22.3  -88.7   59.4   25.4  -54.1   56.5  -48.4  -14.1
## [11]   45.3   -3.9  -30.9   29.6   -3.4  -59.1   89.4  -13.6  -23.3   -7.0
## [21]  -10.6   -1.7   55.6  -65.9   18.5    7.4    3.2   39.9  -37.2    1.0
## [31]   68.1 -119.5   83.3  -83.5  106.6  -20.4  -80.5  116.9  -89.2   33.4
## [41]   -9.5   17.5  -55.3   55.5   17.4  -65.7    9.2   48.6  -43.4   33.4
## [51]   -4.1    2.4   -4.8  -33.0   13.4   70.7  -82.0   29.7

Graficando

Zt1=Zt[3:60]
d1a=d1[2:59]
serdiff=data.frame(Zt1,d1a, d2)
stSeries=ts(serdiff,freq=1,start=1960) #freq. Nuetra serie no tiene estacionalidad
plot(stSeries)

- No se observa una conversion clara de los grados del polinomio debido al componente aleatorio

Pruebas de tendencia

• Comprobar la significancia del componente de tendencia
• Puede ser hecha antes y despues de estimar Tt.
• No deberia estar presente el componente estacional

1. anova

# Ho: La tendencia no es significativa
# H1: La tendencia es significativa
anova(regre2)

## Analysis of Variance Table
## 
## Response: Zt
##           Df Sum Sq Mean Sq  F value Pr(>F)    
## t          1     20      20   0.0323  0.858    
## t2         1 191166  191166 308.2485 <2e-16 ***
## Residuals 57  35350     620                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pruebas no parametricas: para tendencia lineal o o sistematico (no aleatorio)

#cargando datos
serie2=read_excel('serie_tendencia_lineal.xlsx') #tendencia no cuadratica (ascendente o descendente) o sistematico con consecutivo cambio de racha
str(serie2)

## tibble [64 × 2] (S3: tbl_df/tbl/data.frame)
##  $ t : num [1:64] 1 2 3 4 5 6 7 8 9 10 ...
##  $ Zt: num [1:64] 304 294 306 310 304 ...

Z1t=serie2$Zt

2. RACHAS, WALD-WOLFOWITZ

• H0: La secuencia es aleatoria (no hay tendencia)
• H1: La secuencia no es aleatoria (hay tendencia)
• se evalua los datos sobre la mediana, mayor igual le ponemos el simbolo a y menor a la media b
• n1= numeros de signos a y el n2 numero simbolos de n2
• La tabla funciona hasta con 40 datos , con n1 y n2 menor que 20
• suponiendo n1=15, n2=10
• se busca en la tabla 15 y 10 =7 y el otro valor 15 y 10= 18 (valor critico superiores e inferiores)
• Si t <7 pero t >18, se rechaza Ho (t es el estadistico de prueba)
• Cuando n1 o n2>20 se hace una aproximancion a la normal

library(randtests)
runs.test(Z1t,plot = TRUE)

## 
##  Runs Test
## 
## data:  Z1t
## statistic = -2.2681, runs = 24, n1 = 32, n2 = 32, n = 64, p-value =
## 0.02332
## alternative hypothesis: nonrandomness

• hay runs = 24, entonces 26 rachas
• con p-value = 0.02332, se rechaza la hipotesis nula. entonces si hay tendencia
• en el grafico se observan 24 rachas (secuencias)

2. Prueba de signos (COX STUART)

• Ho: La tendencia no es significativa
• H1: La tendencia es significativa
• Agrupa las observaciones en pares
• se puede aproximar una distribucion binomial (para n>20) a la normal con media =n/2, sd=n/4. Con n igual al numero de pares
• Se elimina cuando son iguales

# Z1t: es la serie de los datos
size_=length(Z1t); size_

## [1] 64

Semsize_=size_/2 #cantidad de pares

N=64, entonces c=32 realizando los calculos de camparacion

vec1=Z1t[c(1:Semsize_)]
vec2=Z1t[c(Semsize_+1:Semsize_)]
vec3=vec2-vec1
vec3[vec3>=0] #Si hay cero tener cuidado y eliminar. Se observa que la cantidad de elementos positivos es 4

## [1]  5.0  5.8  1.2 26.1

CantsignosPost=4 #Por lo tanto t2=4, tiene distribucion binomial con n=32 (cantidas de pares, si hubiera algun empate serie menor el n), p=1/2
#En lugar de 4 tambien se podria usar la cantidad de signos negativos
# se le resta cero ya que no hay ninguna igualdad, caso contrario seria 32-1 o 2 o 3 dependiende la cantidad de empates
binom.test(x= CantsignosPost,n=Semsize_ - 0,p=0.5,alternative = "two.sided")

## 
##  Exact binomial test
## 
## data:  CantsignosPost and Semsize_ - 0
## number of successes = 4, number of trials = 32, p-value = 1.93e-05
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.03513065 0.28994842
## sample estimates:
## probability of success 
##                  0.125

• con un p-value = 1.93e-05, se rechaza la hipotesis nula. entonces si hay tendencia

3. Prueba basada en el coeficiente de determinacion de Spearman

• Ordernar la data en funcion a Zt de menor a mayor
• darle rango (1,2,3…)
• verificar empates haciendo Zt - el anterior. Si da cero, hay empate, y se reemplaza el promedio de los dos rangos a cada elemento.
• Luego reordenar por t -hacer la diferencia en otra columna rango - t Luego suma de cuadrados de la nueva comlumna que almacena la diferencia. Este valor es llamado T3

#Zt[1:30]
Rt=rank(serie2$Zt)
Rt

##  [1] 59.0 44.0 61.0 63.0 60.0 57.0 35.0 45.0 55.0  9.0 54.0 51.0 49.0 38.0 27.5
## [16] 47.5 32.0 12.0 21.5 40.5 56.0 62.0 21.5 43.0 33.0 31.0 64.0  7.0 24.0 34.0
## [31] 42.0 40.5  1.0 17.0 46.0 50.0 36.0 52.0  2.0 37.0  3.0 18.0 30.0 47.5 58.0
## [46] 10.5  5.0 20.0 26.0  6.0 25.0 16.0 39.0 29.0 14.0 13.0 27.5  8.0 19.0 53.0
## [61] 23.0 15.0  4.0 10.5

cor.test(serie2$t,Rt,alternative = "two.sided",method = "spearman")

## Warning in cor.test.default(serie2$t, Rt, alternative = "two.sided", method =
## "spearman"): Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  serie2$t and Rt
## S = 67459, p-value = 3.318e-06
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.5443993

• S = 67459 es el T3
• el rho : -0.5443993
• p-value = 3.318e-06. Se rechaza la hipotesis nula, la tendencia es significativa

Estacionalidad

• Cuanto hay una estacionalidad deterministica, cuando es constante, podemos utilizar otra vez el metodo de la regresion pero agregando el componte estacional mediante un calculo mensual de alfa
• Para el calculo del valor estimado no se concidera el termino aleatorio. Por eso en un modelo aditivo podemos sumar el estimado de la tendencia mas el estimado de la estacionalidad (solamente los dos)
• la estacionalidad se repetira anualmente (cada mes)

El modelo cuadratico (regresion):

Estimacion de la tendencia y la estacionalidad (usemos la matriz desarrollada)

serie3=read_excel('MatrizParaEnviarloaR.xlsx')
str(serie3)

## tibble [96 × 14] (S3: tbl_df/tbl/data.frame)
##  $ t  : num [1:96] 1 2 3 4 5 6 7 8 9 10 ...
##  $ t2 : num [1:96] 1 4 9 16 25 36 49 64 81 100 ...
##  $ d1 : num [1:96] 1 0 0 0 0 0 0 0 0 0 ...
##  $ d2 : num [1:96] 0 1 0 0 0 0 0 0 0 0 ...
##  $ d3 : num [1:96] 0 0 1 0 0 0 0 0 0 0 ...
##  $ d4 : num [1:96] 0 0 0 1 0 0 0 0 0 0 ...
##  $ d5 : num [1:96] 0 0 0 0 1 0 0 0 0 0 ...
##  $ d6 : num [1:96] 0 0 0 0 0 1 0 0 0 0 ...
##  $ d7 : num [1:96] 0 0 0 0 0 0 1 0 0 0 ...
##  $ d8 : num [1:96] 0 0 0 0 0 0 0 1 0 0 ...
##  $ d9 : num [1:96] 0 0 0 0 0 0 0 0 1 0 ...
##  $ d10: num [1:96] 0 0 0 0 0 0 0 0 0 1 ...
##  $ d11: num [1:96] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Zt : num [1:96] 480 469 466 449 450 ...

Regresion lineal para estimar la tendencia y estacionalidad

regre2s=lm(data =serie3,Zt~.)
regre2s

## 
## Call:
## lm(formula = Zt ~ ., data = serie3)
## 
## Coefficients:
## (Intercept)            t           t2           d1           d2           d3  
##   448.06009      5.19641     -0.06289     41.15705     21.20696      7.83265  
##          d4           d5           d6           d7           d8           d9  
##   -18.56588    -29.56363    -37.06060    -36.54430    -18.35222      0.29064  
##         d10          d11  
##    -4.27823     30.76619

summary(regre2s)

## 
## Call:
## lm(formula = Zt ~ ., data = serie3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -46.436 -13.971  -0.847  14.020  64.159 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 448.060089   7.189889  62.318  < 2e-16 ***
## t             5.196413   0.341442  15.219  < 2e-16 ***
## t2           -0.062889   0.003409 -18.449  < 2e-16 ***
## d1           41.157046   7.777595   5.292 9.90e-07 ***
## d2           21.206960   7.772735   2.728  0.00778 ** 
## d3            7.832653   7.768953   1.008  0.31632    
## d4          -18.565877   7.766178  -2.391  0.01911 *  
## d5          -29.563630   7.764357  -3.808  0.00027 ***
## d6          -37.060604   7.763456  -4.774 7.79e-06 ***
## d7          -36.544301   7.763456  -4.707 1.01e-05 ***
## d8          -18.352221   7.764357  -2.364  0.02046 *  
## d9            0.290637   7.766178   0.037  0.97024    
## d10          -4.278227   7.768953  -0.551  0.58335    
## d11          30.766187   7.772735   3.958  0.00016 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.93 on 82 degrees of freedom
## Multiple R-squared:  0.8776, Adjusted R-squared:  0.8581 
## F-statistic: 45.21 on 13 and 82 DF,  p-value: < 2.2e-16

coefs=as.data.frame(summary(regre2s)$coefficients)
coefs$Estimate=round(coefs$Estimate,3)#redondeando a 3 decimales
coefs$Estimate

##  [1] 448.060   5.196  -0.063  41.157  21.207   7.833 -18.566 -29.564 -37.061
## [10] -36.544 -18.352   0.291  -4.278  30.766

• Los tres primeros coeficientes son los que utilizaran para realizar el calculo de la tendencia: 448.060 5.196 -0.063
• Con suma de los siguientes coeficientes (41.157 21.207 7.833 -18.566 -29.564 -37.061 -36.544 -18.352 0.291 -4.278 30.766) se obtiene el d12

Obteniendo coeficiente 12

# solo los d's
vectorEfec=coefs$Estimate
sumar_=vectorEfec[c(4:length(vectorEfec))]
sumar_

##  [1]  41.157  21.207   7.833 -18.566 -29.564 -37.061 -36.544 -18.352   0.291
## [10]  -4.278  30.766

d12=sum(sumar_)*-1
d12

## [1] 43.111

min_xlim=1
max_xlim=96
#creando los d's
numero = 0
d1 = rep(numero, times = max_xlim)
d2=d1
d3=d1
d4=d1
d5=d1
d6=d1
d7=d1
d8=d1
d9=d1
d10=d1
d11=d1
#Prediciendo la tendencia
datan=data.frame(t=seq(min_xlim,max_xlim),
                 t2=(seq(min_xlim,max_xlim)**2))
datan$Zt=serie3$Zt #valor real
datan$Tt=predict(regre2s,newdata = datan)

#los d's (todos los d's) son 12
ds=sumar_
ds = c(ds, d12)
ds

##  [1]  41.157  21.207   7.833 -18.566 -29.564 -37.061 -36.544 -18.352   0.291
## [10]  -4.278  30.766  43.111

#agregaremos al df el componente de tendencia
repeticiones <- 8
vector_repetido <- rep(ds, times = repeticiones)
datan$St=vector_repetido

#Adicionado el componente aleatorio
datan$at=datan$Zt - (datan$Tt+datan$St)
#calculado el valor de la serie estimado.
datan$Ztest=datan$Tt+datan$St
head(datan,5)

##   t t2    Zt       Tt      St          at    Ztest
## 1 1  1 480.4 453.1936  41.157 -13.9506128 494.3506
## 2 2  4 468.7 458.2014  21.207 -10.7083588 479.4084
## 3 3  9 465.5 463.0833   7.833  -5.4163273 470.9163
## 4 4 16 448.7 467.8395 -18.566  -0.5735181 449.2735
## 5 5 25 449.8 472.4699 -29.564   6.8940688 442.9059

Significancia de la estacionalidad

• H0:α1=α2=…=α12=0 estacionalidad no significativa
• H1: αi≠0, algún i, estacionalidad significativa

regre2s_reducido = lm(data =serie3,Zt~serie3$t+serie3$t2) #modelo reducido
anova(regre2s_reducido)

## Analysis of Variance Table
## 
## Response: Zt
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## serie3$t   1  59131   59131   46.82 8.145e-10 ***
## serie3$t2  1 175649  175649  139.08 < 2.2e-16 ***
## Residuals 93 117453    1263                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(regre2s)

## Analysis of Variance Table
## 
## Response: Zt
##           Df Sum Sq Mean Sq  F value    Pr(>F)    
## t          1  59131   59131 112.4214 < 2.2e-16 ***
## t2         1 175649  175649 333.9500 < 2.2e-16 ***
## d1         1     14      14   0.0266 0.8708432    
## d2         1   2331    2331   4.4313 0.0383464 *  
## d3         1   4479    4479   8.5154 0.0045433 ** 
## d4         1  14070   14070  26.7505 1.608e-06 ***
## d5         1  15684   15684  29.8181 4.953e-07 ***
## d6         1  15746   15746  29.9367 4.736e-07 ***
## d7         1  11525   11525  21.9122 1.115e-05 ***
## d8         1   2139    2139   4.0659 0.0470296 *  
## d9         1     74      74   0.1407 0.7085591    
## d10        1     21      21   0.0399 0.8421139    
## d11        1   8241    8241  15.6675 0.0001599 ***
## Residuals 82  43130     526                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• la suma de cuadrados de error del modelo reducido es: 117453, con 93 GL
• la suma de cuadrados de error del modelo completo es: 43130, con 82 GL

Contribucion del componente estacional (haciendo la diferencia)

scH= 117453- 43130
dGL=93-82
scH

## [1] 74323

dGL

## [1] 11

• La contribucion del componente estacional(H) es: 74323, GL= 11

El F calculado (Fc)

Fc = (74323/11)/(43130/82) #El denominador son los datos del modelo completo -- 11, 82
Fc

## [1] 12.84591

Obteniendo el P-value

pf(12.84591,11,82, lower.tail = F)

## [1] 9.106623e-14

• con p-value de 9.106623e-14, se rechaza la hipotesis nula. si hay significancia de la tendencia

El metodo de medias moviles

• recordemos que decomponse usa medias moviles

serie4=read_excel('serie_anio_mes_zt.xlsx')
str(serie4)

## tibble [96 × 5] (S3: tbl_df/tbl/data.frame)
##  $ t  : num [1:96] 1 2 3 4 5 6 7 8 9 10 ...
##  $ t2 : num [1:96] 1 4 9 16 25 36 49 64 81 100 ...
##  $ MES: num [1:96] 1 2 3 4 5 6 7 8 9 10 ...
##  $ AÑO: num [1:96] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Zt : num [1:96] 480 469 466 449 450 ...

Usando la funcion decompose: Primero dar el formato de serie

Zt=serie4$Zt
Zts=ts(Zt,frequency = 12, start=c(2012,1)) #nuestra serie esta en meses (12)
decom=decompose(Zts,type = 'additive') #para un modelo multiplicative
#decom
#$x son los datos de la serie
#$figure indices estacionales (inicia en enero)

Tendencia

decom$trend

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2012       NA       NA       NA       NA       NA       NA 478.6458 482.4333
## 2013 499.6625 503.9792 508.5333 514.3917 518.4792 522.7833 527.0458 530.8125
## 2014 544.1958 544.4458 543.1583 543.5333 548.1583 551.3250 552.8792 555.6375
## 2015 557.5875 559.7000 560.7792 556.6792 550.2125 547.7708 547.7292 545.9583
## 2016 546.3333 543.7208 542.3833 543.1125 545.6042 546.9208 546.6250 545.8333
## 2017 533.9292 529.1250 524.6167 522.0083 519.8875 514.9792 509.4125 505.4750
## 2018 489.2583 485.1708 480.4375 473.9500 467.8333 463.1917 459.6708 456.8833
## 2019 431.8792 428.7583 425.1542 421.9292 415.7750 410.7583       NA       NA
##           Sep      Oct      Nov      Dec
## 2012 485.4250 489.5417 492.0125 494.3583
## 2013 534.4292 535.5875 539.1542 543.5833
## 2014 557.6458 559.6875 560.4125 558.4917
## 2015 546.7583 546.0583 544.5250 546.4125
## 2016 542.1625 541.6458 540.6667 536.8042
## 2017 502.8792 498.4125 495.5125 493.9417
## 2018 453.6208 450.4458 444.8875 436.4917
## 2019       NA       NA       NA       NA

Indices estacionales

decom$figure

##  [1]  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
##  [7] -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417

Componente estacional

decom$seasonal

##              Jan         Feb         Mar         Apr         May         Jun
## 2012  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2013  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2014  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2015  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2016  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2017  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2018  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
## 2019  43.0421131  22.7343750   8.8540179 -18.1662202 -30.2584821 -37.2840774
##              Jul         Aug         Sep         Oct         Nov         Dec
## 2012 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2013 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2014 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2015 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2016 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2017 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2018 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417
## 2019 -37.8525298 -16.6846726  -0.4971726  -4.3769345  32.0385417  38.4510417

Componente aleatorio

decom$random

##                Jan           Feb           Mar           Apr           May
## 2012            NA            NA            NA            NA            NA
## 2013 -1.204613e+00 -2.821354e+01 -9.887351e+00  9.274554e+00 -3.592068e+01
## 2014 -1.003795e+01 -1.398021e+01 -1.241235e+01 -2.416711e+01  2.430015e+01
## 2015  1.270387e+00  1.226563e+01 -2.333318e+01  4.987054e+00 -2.654018e+00
## 2016 -1.197545e+01  1.024479e+01  3.226265e+01 -3.544628e+01  1.915432e+01
## 2017  4.628720e+00  1.640625e+00 -1.487068e+01  3.815789e+01 -3.112902e+01
## 2018 -4.464286e-04  3.947917e-01  1.220848e+01 -3.883780e+00  4.142515e+01
## 2019  2.347872e+01  2.380729e+01  2.219182e+01  1.723705e+01 -9.016518e+00
##                Jun           Jul           Aug           Sep           Oct
## 2012            NA  1.370670e+01  1.375134e+01  2.977217e+01 -6.164732e+00
## 2013  9.300744e+00  3.880670e+01 -4.527827e+00  5.996801e+01  9.189435e+00
## 2014 -2.840923e+00  1.127336e+01 -2.165283e+01 -1.848661e+00  3.268943e+01
## 2015 -2.048676e+01  1.592336e+01  3.922634e+01 -1.626116e+01 -2.678140e+01
## 2016  8.463244e+00 -1.297247e+01  6.651339e+00 -1.106533e+01 -5.468899e+00
## 2017  2.370491e+01 -2.805997e+01 -1.599033e+01 -1.698199e+01  2.036443e+01
## 2018  1.729241e+01 -3.251830e+01 -1.129866e+01 -3.742366e+01 -1.766890e+01
## 2019 -2.927426e+01            NA            NA            NA            NA
##                Nov           Dec
## 2012 -4.351042e+00 -2.110938e+01
## 2013 -1.479271e+01 -3.734375e+00
## 2014  2.734896e+01 -6.042708e+00
## 2015 -3.886354e+01  2.953646e+01
## 2016  7.894792e+00  2.784479e+01
## 2017  1.954896e+01 -1.359271e+01
## 2018  9.373958e+00 -6.742708e+00
## 2019            NA            NA

Juntando en un solo dataframe todo

serie4$Tt=as.vector(decom$trend)
serie4$St=as.vector(decom$seasonal)
serie4$at=as.vector(decom$random)
serie4$Ztest=serie4$Tt+serie4$St #para un modelo aditivo

head(serie4,10)

## # A tibble: 10 × 9
##        t    t2   MES   AÑO    Zt    Tt      St    at Ztest
##    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>
##  1     1     1     1     1  480.   NA   43.0   NA      NA 
##  2     2     4     2     1  469.   NA   22.7   NA      NA 
##  3     3     9     3     1  466.   NA    8.85  NA      NA 
##  4     4    16     4     1  449.   NA  -18.2   NA      NA 
##  5     5    25     5     1  450.   NA  -30.3   NA      NA 
##  6     6    36     6     1  441    NA  -37.3   NA      NA 
##  7     7    49     7     1  454.  479. -37.9   13.7   441.
##  8     8    64     8     1  480.  482. -16.7   13.8   466.
##  9     9    81     9     1  515.  485.  -0.497 29.8   485.
## 10    10   100    10     1  479   490.  -4.38  -6.16  485.

• Guardando en excel los resultados

library("writexl")
write_xlsx(serie4,"output.xlsx")  

Para un modelo multiplicativo Para un modelo multiplicativo el calculo de la tendencia es la misma

Prueba de significancia del componente estacional por

La prueba de KRUSKALL WALLIS (prueba no parametrica)

• Esta pruaba es no parametrica
• Ho: La tendencia no es significativa
• H1: La tendencia es significativa
• Esta prueba trabaja ranqueando Yt

library(car)

## Loading required package: carData

año=serie4$AÑO
mes=serie4$MES
Z3t=serie4$Zt
Yt=Z3t-decom$trend # decom componente aditivo creado antes 
#Yt
kruskal.test(Yt ~ mes, data = serie4)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Yt by mes
## Kruskal-Wallis chi-squared = 54.134, df = 11, p-value = 1.116e-07

• con un p-value 1.116e-07 se rechaza la hipotesis nula, . En conclusion si hay tendencia

Prueba de friedman WALLIS (prueba no parametrica)

• Esta pruaba es no parametrica
• Ho: La tendencia no es significativa
• H1: La tendencia es significativa
• Esta prueba trabaja ranqueando Zt

añof=as.factor(año)
mesf=as.factor(mes)
friedman.test(Z3t,mesf,añof)

## 
##  Friedman rank sum test
## 
## data:  Z3t, mesf and añof
## Friedman chi-squared = 47.577, df = 11, p-value = 1.698e-06

• con un p-value 1.698e-06 se rechaza la hipotesis nula, . En conclusion si hay tendencia

prueba de significancia por Anova (parametricas)

#1 VÍA: se usa Yt
anova(lm(Yt ~ mesf))

## Analysis of Variance Table
## 
## Response: Yt
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## mesf      11  65234  5930.3  11.333 7.206e-12 ***
## Residuals 72  37677   523.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• se observa significancia del mes, lo cual indica existencia de estacionalidad

#2 VÍA: se usa Zt
anova(lm(Z3t ~ mesf+añof))

## Analysis of Variance Table
## 
## Response: Z3t
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## mesf      11  70802  6436.5  7.7717 7.135e-09 ***
## añof       7 217660 31094.2 37.5446 < 2.2e-16 ***
## Residuals 77  63771   828.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• se observa significancia del mes, lo cual indica existencia de estacionalidad (aqui importa el mes)