Simulaciones¶

Autorregresivo de orden 1 y 2¶

$$ Z_{t}\sim Normal(\mu, \sigma^{2}) $$
In [4]:
# Ruido blanco para t=0,...,100
# mu = 0
# sigma^2 = 1--->sigma(desviacion estandar) = 1
Zt <- rnorm(100, 0, 1)

ts.plot(Zt, main="Ruido blanco")
grid()

Autorregresivo de orden 1¶

$$ X_{t}=\alpha_{0}+\alpha_{1}X_{t-1}+Z_{t} $$
$$ X_{1}=2-0.8X_{0}+Z_{1} $$
$$ X_{1}=2-0.8(0)-1.2424=0.75 $$
$$ X_{2}=2-0.8(0.75)-0.44 $$
In [10]:
Zt[2]
-0.442100553519466
In [11]:
v1 <- c(0)
# acceder al primer elemento del vector
v1[2] <- 7
v1[3] <- 700
v1[4] <- 7000

v1
  1. 0
  2. 7
  3. 700
  4. 7000
In [17]:
# Alpha
a0 = 2
# alpha1 --> (-1, 1)
a1 = 0.5

# inicializar el modelo
# t=0
Xt <- c(0)

# rellenar el resto de valores Xt

for (t in 2:100){
    Xt[t] = a0 + a1 * Xt[t-1] + Zt[t]
}

ts.plot(Xt, main="AR(1)")
$$ X_{t}=\alpha_{0}+\alpha_{1}X_{t-1}+Z_{t} $$
In [18]:
# Alpha
a0 = 2
# alpha1 --> (-1, 1)
a1 = 0.5
a2 = -0.8

# inicializar el modelo
# t=0
Xt <- c(0,0)

# rellenar el resto de valores Xt

for (t in 3:100){
    Xt[t] = a0 + a1 * Xt[t-1] + a2 * Xt[t-2] + Zt[t]
}

ts.plot(Xt, main="AR(2)")
$$ X_{t}=Z_{t}+\beta_{1}Z_{t-1} $$
In [19]:
# beta1 (-1,1)
b1 = 0.5

# inicializar el modelo
# t=0
Xt <- c(0)

# rellenar el resto de valores Xt

for (t in 2:100){
    Xt[t] = Zt[t] + b1 * Zt[t-1]
}

ts.plot(Xt, main="MA(1)")
$$ X_{t}=\alpha_{0}+\alpha_{1}X_{t-1}+Z_{t}+\beta_{1}Z_{t-1} $$
In [20]:
# alpha's
a0 = 2
a1 = -0.8

# beta1 (-1,1)
b1 = 0.5

# inicializar el modelo
# t=0
Xt <- c(0)

# rellenar el resto de valores Xt

for (t in 2:100){
    Xt[t] = a0 + a1 * Xt[t-1] + Zt[t] + b1 * Zt[t-1]
}

ts.plot(Xt, main="ARMA(p=1,q=1)")

ARIMA(p,d,q):¶

In [21]:
#                             p,d,q   alpha1   beta1
xt <- arima.sim(model= list(c(1,0,1), ar=-0.9, ma=0.5) , n=100)
ts.plot(xt, main="ARMA(1,1), alpha1=-0.9, beta1=0.5")
In [23]:
# Si vemos tendencia y ciclo estacional:
# d=2
#                             p,d,q   alpha1   beta1
#xt <- arima.sim(model= list(c(1,2,1), ar=-0.9, ma=0.5) , n=100)
#ts.plot(xt, main="ARMA(1,1), alpha1=-0.9, beta1=0.5")
In [24]:
# Consideremos
Xt <- AirPassengers
Xt
A Time Series: 12 × 12
JanFebMarAprMayJunJulAugSepOctNovDec
1949112118132129121135148148136119104118
1950115126141135125149170170158133114140
1951145150178163172178199199184162146166
1952171180193181183218230242209191172194
1953196196236235229243264272237211180201
1954204188235227234264302293259229203229
1955242233267269270315364347312274237278
1956284277317313318374413405355306271306
1957315301356348355422465467404347305336
1958340318362348363435491505404359310337
1959360342406396420472548559463407362405
1960417391419461472535622606508461390432
In [25]:
ts.plot(Xt)
In [27]:
# importacion necesaria
library(forecast)

fit <- auto.arima(Xt)
fit

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 


Series: Xt 
ARIMA(2,1,1)(0,1,0)[12] 

Coefficients:
         ar1     ar2      ma1
      0.5960  0.2143  -0.9819
s.e.  0.0888  0.0880   0.0292

sigma^2 = 132.3:  log likelihood = -504.92
AIC=1017.85   AICc=1018.17   BIC=1029.35
In [28]:
# pronostico a 12 tiempos futuros (un año)
pronostico <- forecast(fit,12,level=95)

# graficamos
plot(pronostico, main="Pronóstico con auto.arima()")
grid()
In [29]:
pronosticos_t <- data.frame(pronostico$mean,pronostico$lower,pronostico$upper)
pronosticos_t
A data.frame: 12 × 3
pronostico.meanX95.X95..1
<ts><dbl><dbl>
445.6349423.0851468.1847
420.3950393.9304446.8596
449.1983419.4892478.9074
491.8399460.0092523.6707
503.3945469.9953536.7937
566.8625532.3007601.4242
654.2602618.8122689.7081
638.5975602.4630674.7320
540.8837504.2081577.5594
494.1266457.0177531.2356
423.3327385.8715460.7940
465.5076427.7556503.2596
In [30]:
# Xt: Ts, St, Zt
plot(decompose(Xt))