Моделирование данных временных рядов с пояснительной переменной - PullRequest
1 голос
/ 20 июня 2019

У меня есть данные ниже, где

  • у - средняя цена корзины некоторых продуктов
  • х - это обменный курс.
   period    y    x
1  201501 1530 2.49
2  201502 1450 2.62
3  201503 1637 2.77
4  201504 1404 2.84
5  201505 1442 2.82
6  201506 1442 2.89
7  201507 1518 2.88
8  201508 1492 3.05
9  201509 1743 3.21
10 201510 1902 3.14
11 201511 1855 3.07
12 201512 1879 3.12
13 201601 2018 3.21
14 201602 2117 3.15
15 201603 2002 3.09
16 201604 1837 3.04
17 201605 1902 3.14
18 201606 1910 3.12
19 201607 2162 3.16
20 201608 2183 3.17
21 201609 2100 3.17
22 201610 2122 3.28
23 201611 2461 3.51
24 201612 2250 3.73
25 201701 2466 4.00
26 201702 2212 3.93
27 201703 2424 3.93
28 201704 2477 3.91
29 201705 2402 3.82
30 201706 2360 3.77
31 201707 2475 3.81
32 201708 2690 3.76
33 201709 2655 3.70
34 201710 2889 3.92
35 201711 2683 4.15
36 201712 2674 4.12
37 201801 2695 4.03
38 201802 2707 4.04
39 201803 2728 4.15
40 201804 2607 4.33
41 201805 2917 4.71
42 201806 2946 4.94
43 201807 3031 5.08
44 201808 3224 6.20
45 201809 3962 6.76
46 201810 4043 6.25
47 201811 3805 5.76
48 201812 3607 5.67
49 201901 3694 5.74
50 201902 3566 5.63
51 201903 3541 5.83
52 201904 3350 6.15
  • y, как полагают, сильно зависит от x и его лагов
  • у вас есть подозрения на сезонность.
  • возможно, на y влияют другие переменные, но эти переменные недоступны.

Итак, чтобы выявить связь между у тревогой х (и отставанием от х), я сделал следующий анализ:

library(DataCombine)

#data
data2<-structure(c(201501, 201502, 201503, 201504, 201505, 201506, 201507,
                 201508, 201509, 201510, 201511, 201512, 201601, 201602, 201603,
                 201604, 201605, 201606, 201607, 201608, 201609, 201610, 201611,
                 201612, 201701, 201702, 201703, 201704, 201705, 201706, 201707,
                 201708, 201709, 201710, 201711, 201712, 201801, 201802, 201803,
                 201804, 201805, 201806, 201807, 201808, 201809, 201810, 201811,
                 201812, 201901, 201902, 201903, 201904, 1530, 1450, 1637, 1404,
                 1442, 1442, 1518, 1492, 1743, 1902, 1855, 1879, 2018, 2117, 2002,
                 1837, 1902, 1910, 2162, 2183, 2100, 2122, 2461, 2250, 2466, 2212,
                 2424, 2477, 2402, 2360, 2475, 2690, 2655, 2889, 2683, 2674, 2695,
                 2707, 2728, 2607, 2917, 2946, 3031, 3224, 3962, 4043, 3805, 3607,
                 3694, 3566, 3541, 3350, 2.49, 2.62, 2.77, 2.84, 2.82, 2.89, 2.88,
                 3.05, 3.21, 3.14, 3.07, 3.12, 3.21, 3.15, 3.09, 3.04, 3.14, 3.12,
                 3.16, 3.17, 3.17, 3.28, 3.51, 3.73, 4, 3.93, 3.93, 3.91, 3.82,
                 3.77, 3.81, 3.76, 3.7, 3.92, 4.15, 4.12, 4.03, 4.04, 4.15, 4.33,
                 4.71, 4.94, 5.08, 6.2, 6.76, 6.25, 5.76, 5.67, 5.74, 5.63, 5.83,
                 6.15), .Dim = c(52L, 3L), .Dimnames = list(NULL, c("period", 
                 "y", "x")), .Tsp = c(1, 5.25, 12), 
                 class = c("mts", "ts", "matrix"))

# deaseasonal data using loess procedure
# model assumed to be multiplicative ->
# so seasonal coefficients obtained after taking logarithm


data2<-ts(data1, frequency = 12)
lprod<-log(data2[,2])
decomp<-stl(lprod, s.window="periodic")

decomp<-decomp$time.series

season<-exp(decomp[,1])
trend<-exp(decomp[,2])
rand<-exp(decomp[,3])

#deasonal y value
desdata<-trend*rand


#obtaining lags(1 and 2) of explonatary variable x
ex_var<-as.data.frame(data2)
ex_var<-slide(ex_var, Var='x', NewVar = "x1", slideBy = -1)  
ex_var<-slide(ex_var, "x", NewVar = "x2", slideBy = -2)  
ex_var<-slide(ex_var, "x", NewVar = "x3", slideBy = -2)  

#delete firts two rows
ex_var<-ex_var[-c(1:2),]

#regression
#I also include x variable at time t. My aim is just to obtaion the relation
#After some trials, the below models is fitted

myreg<-lm(formula=y~-1+x+x1, data=ex_var)

#fitted values and deviations
fitted<-myreg$fitted.values
dev<-(fitted/ex_var$y-1)*100

Итак, я просто десезонализировал y и регрессировал y на x и x1.

Затем я получил установленные значения и отклонения.

     dev
3    1.8692282
4   24.5705019
5   23.3397058
6   23.4631450
7   19.3359159
8   22.9155774
9   11.2428038
10   5.2958604
11   5.5962845
12   2.9269927
13  -2.2926331
14  -5.3226638
15  -1.7643976
16   5.0966933
17   1.1112870
18   2.9722925
19  -9.1685959
20  -9.1115852
21  -5.2963860
22  -5.4530041
23 -14.8963143
24  -0.5730753
25  -3.3622134
26  12.9454251
27   1.7153551
28  -0.5895710
29   1.5281120
30   1.2122606
31  -4.1790659
32 -11.4373380
33 -11.5113009
34 -18.4386679
35  -6.9727496
36  -2.8112683
37  -4.6213799
38  -6.5419446
39  -6.4478475
40   0.9688325
41  -4.7986412
42   1.5460511
43   2.9863159
44   4.3845503
45   0.4257484
46   2.8904249
47   1.0008665
48  -0.2121615
49  -3.4013337
50   0.4939980
51   0.6482636
52  10.7024816

Кажется, что отклонения сильно автокоррелированы.

Итак, я думаю, мне нужно построить модель, которая включает в себя отставание от y. Но я не уверен в этом. Нужно ли принимать различия логарифмов? Или я должен идти другим путем?

Я не эксперт по данным временных рядов. Я действительно застрял с этими данными. Я буду очень рад за любую помощь. Большое спасибо за то, что вы уделили мне время.

1 Ответ

1 голос
/ 21 июня 2019

Я не могу определить сезонность.
Для линейной регрессии с лагами вы можете использовать dynlm, но вы должны также рассмотреть ( псевдо ) модель ARIMAX.

library(zoo)
library(dynlm)

# Convert to a multivariate zoo object
ym <- as.yearmon(as.character(data2[,"period"]), "%Y%m")
tt.zoo <- zoo(data2[,c("y", "x")], order.by=ym)

# No significant periodicity
tt.d <- diff(tt.zoo)
ccf(tt.d[,"x"], tt.d[,"y"])
acf(tt.d[,"y"])
acf(tt.d[,"x"])

# Dynamic Linear Models of various complexity
dlm0 <- dynlm(y ~ L(y, 1), 
  data=tt.zoo, start=start(tt.zoo)+1/12, end=end(tt.zoo))
dlm1 <- dynlm(y ~ x, 
  data=tt.zoo, start=start(tt.zoo)+1/12, end=end(tt.zoo))
dlm2 <- dynlm(y ~ x + L(x, 1), 
  data=tt.zoo, start=start(tt.zoo)+1/12, end=end(tt.zoo))
dlm3 <- dynlm(y ~ x + L(x, 1) + L(y, 1), 
  data=tt.zoo, start=start(tt.zoo)+1/12, end=end(tt.zoo))

# Residuals look reasonably OK. Maybe a slight curve?
plot(residuals(dlm0), col=2)
lines(residuals(dlm1), col=3)
lines(residuals(dlm2), col=4)
lines(residuals(dlm3), col=5)

# Fits are pretty OK
plot(tt.zoo[,"y"], col="grey", lwd=2)
lines(fitted(dlm0), col=2)
lines(fitted(dlm1), col=3)
lines(fitted(dlm2), col=4)
lines(fitted(dlm3), col=5)

# dlm1 and dlm2 have significant acf, due to y's non-stationarity
acf(residuals(dlm0))
acf(residuals(dlm1))
acf(residuals(dlm2))
acf(residuals(dlm3))

# dlm3 seems like it makes good use of the extra marameters
anova(dlm0, dlm1, dlm2, dlm3)
AIC(dlm0, dlm1, dlm2, dlm3)
#      df      AIC
# dlm0  3 678.6580
# dlm1  3 699.1142
# dlm2  4 686.8249
# dlm3  5 660.2720


# Altermative 'ARIMAX' model with lagged external regressor
tt.l <- cbind(tt.zoo, xl=lag(tt.zoo[,2], -1))
tt.l <- tt.l[complete.cases(tt.l),]

ax1 <- arima(tt.l[,1], order=c(0, 1, 0), xreg=tt.l[,2])
ax2 <- arima(tt.l[,1], order=c(0, 1, 0), xreg=tt.l[,2:3])
ax3 <- arima(tt.l[,1], order=c(1, 1, 0), xreg=tt.l[,2:3])

# Again the extra parameters seems to be justified
AIC(ax1, ax2, ax3)
#     df      AIC
# ax1  2 656.1224
# ax2  3 649.2153
# ax3  4 645.8715

# Checking that the residuals are sufficiently white noise-like
library(forecast)
checkresiduals(ax1)
checkresiduals(ax2)
checkresiduals(ax3)
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...