1) Jeśli chodzi o twoje pierwsze pytanie, w literaturze opracowano i omówiono niektóre statystyki testów w celu przetestowania zerowej stacjonarności i zerowej wartości pierwiastka jednostkowego. Niektóre z wielu artykułów napisanych na ten temat to:
Związane z trendem:
- Dickey, D. y Fuller, W. (1979a), Rozkład estymatorów dla autoregresywnych szeregów czasowych z pierwiastkiem jednostkowym, Journal of American Statistics Association 74, 427-31.
- Dickey, D. y Fuller, W. (1981), Statystyka ilorazu wiarygodności dla szeregów czasowych autoregresji z pierwiastkiem jednostkowym, Econometrica 49, 1057-1071.
- Kwiatkowski, D., Phillips, P., Schmidt, P. y Shin, Y. (1992), Testowanie zerowej hipotezy stacjonarności wobec alternatywy pierwiastka: Czy jesteśmy pewni, że ekonomiczne szeregi czasowe mają pierwiastek? , Journal of Econometrics 54, 159-178.
- Phillips, P. y Perron, P. (1988), Testowanie pierwiastka jednostkowego w regresji szeregów czasowych, Biometrika 75, 335-46.
- Durlauf, S. y Phillips, P. (1988), Trendy a losowe spacery w analizie szeregów czasowych, Econometrica 56, 1333-54.
Związane ze składnikiem sezonowym:
- Hylleberg, S., Engle, R., Granger, C. y Yoo, B. (1990), Sezonowa integracja i kointegracja, Journal of Econometrics 44, 215-38.
- Canova, F. y Hansen, BE (1995), Czy wzorce sezonowe są stałe w czasie? test stabilności sezonowej, Journal of Business and Economic Statistics 13, 237-252.
- Franses, P. (1990), Testowanie sezonowych korzeni jednostkowych w danych miesięcznych, Raport techniczny 9032, Instytut Ekonometryczny.
- Ghysels, E., Lee, H. y Noh, J. (1994), Testowanie pierwiastków jednostkowych w sezonowych szeregach czasowych. niektóre rozszerzenia teoretyczne i dochodzenie Monte Carlo, Journal of Econometrics 62, 415-442.
Podręcznik Banerjee, A., Dolado, J., Galbraith, J. y Hendry, D. (1993), Kointegracja, korekcja błędów oraz analiza ekonometryczna danych niestacjonarnych, Advanced Texts in Econometrics. Oxford University Press jest również dobrym punktem odniesienia.
2) Twoja druga troska jest uzasadniona literaturą. Jeśli istnieje test na pierwiastek jednostkowy, tradycyjna statystyka t, którą zastosowałbyś do trendu liniowego, nie jest zgodna z rozkładem standardowym. Patrz na przykład Phillips, P. (1987), Regresja szeregów czasowych z pierwiastkiem jednostkowym, Econometrica 55 (2), 277-301.
Jeśli pierwiastek jednostkowy istnieje i jest ignorowany, wówczas zmniejsza się prawdopodobieństwo odrzucenia wartości zerowej, że współczynnik trendu liniowego wynosi zero. Oznacza to, że zbyt często skończylibyśmy na modelowaniu deterministycznego trendu liniowego dla danego poziomu istotności. W obecności korzenia jednostki powinniśmy zamiast tego przekształcić dane, biorąc regularne różnice do danych.
3) Na przykład, jeśli używasz R, możesz wykonać następującą analizę z danymi.
x <- structure(c(7657, 5451, 10883, 9554, 9519, 10047, 10663, 10864,
11447, 12710, 15169, 16205, 14507, 15400, 16800, 19000, 20198,
18573, 19375, 21032, 23250, 25219, 28549, 29759, 28262, 28506,
33885, 34776, 35347, 34628, 33043, 30214, 31013, 31496, 34115,
33433, 34198, 35863, 37789, 34561, 36434, 34371, 33307, 33295,
36514, 36593, 38311, 42773, 45000, 46000, 42000, 47000, 47500,
48000, 48500, 47000, 48900), .Tsp = c(1, 57, 1), class = "ts")
Najpierw możesz zastosować test Dickeya-Fullera dla wartości zerowej katalogu głównego:
require(tseries)
adf.test(x, alternative = "explosive")
# Augmented Dickey-Fuller Test
# Dickey-Fuller = -2.0685, Lag order = 3, p-value = 0.453
# alternative hypothesis: explosive
oraz test KPSS dla hipotezy odwrotnego zerowania, stacjonarność wobec alternatywy stacjonarności wokół trendu liniowego:
kpss.test(x, null = "Trend", lshort = TRUE)
# KPSS Test for Trend Stationarity
# KPSS Trend = 0.2691, Truncation lag parameter = 1, p-value = 0.01
Wyniki: test ADF, przy poziomie istotności 5% pierwiastek jednostkowy nie jest odrzucany; Test KPSS, zerowa stacjonarność jest odrzucana na rzecz modelu z trendem liniowym.
Na marginesie: użycie lshort=FALSE
wartości zerowej testu KPSS nie jest odrzucane na poziomie 5%, jednak wybiera 5 opóźnień; kolejna kontrola, której nie pokazano, sugeruje, że wybór 1-3 opóźnień jest odpowiedni dla danych i prowadzi do odrzucenia hipotezy zerowej.
Zasadniczo powinniśmy kierować się testem, dla którego byliśmy w stanie odrzucić hipotezę zerową (zamiast testem, dla którego nie odrzuciliśmy (przyjęliśmy) zerową). Jednak regresja oryginalnej serii według trendu liniowego okazuje się niewiarygodna. Z jednej strony kwadrat R jest wysoki (ponad 90%), co jest wskazane w literaturze jako wskaźnik fałszywej regresji.
fit <- lm(x ~ 1 + poly(c(time(x))))
summary(fit)
#Coefficients:
# Estimate Std. Error t value Pr(>|t|)
#(Intercept) 28499.3 381.6 74.69 <2e-16 ***
#poly(c(time(x))) 91387.5 2880.9 31.72 <2e-16 ***
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 2881 on 55 degrees of freedom
#Multiple R-squared: 0.9482, Adjusted R-squared: 0.9472
#F-statistic: 1006 on 1 and 55 DF, p-value: < 2.2e-16
Z drugiej strony reszty są autokorelowane:
acf(residuals(fit)) # not displayed to save space
Co więcej, zerowy rdzeń jednostkowy w resztach nie może zostać odrzucony.
adf.test(residuals(fit))
# Augmented Dickey-Fuller Test
#Dickey-Fuller = -2.0685, Lag order = 3, p-value = 0.547
#alternative hypothesis: stationary
W tym momencie możesz wybrać model, który będzie używany do uzyskiwania prognoz. Na przykład prognozy oparte na strukturalnym modelu szeregów czasowych i na modelu ARIMA można uzyskać w następujący sposób.
# StructTS
fit1 <- StructTS(x, type = "trend")
fit1
#Variances:
# level slope epsilon
#2982955 0 487180
#
# forecasts
p1 <- predict(fit1, 10, main = "Local trend model")
p1$pred
# [1] 49466.53 50150.56 50834.59 51518.62 52202.65 52886.68 53570.70 54254.73
# [9] 54938.76 55622.79
# ARIMA
require(forecast)
fit2 <- auto.arima(x, ic="bic", allowdrift = TRUE)
fit2
#ARIMA(0,1,0) with drift
#Coefficients:
# drift
# 736.4821
#s.e. 267.0055
#sigma^2 estimated as 3992341: log likelihood=-495.54
#AIC=995.09 AICc=995.31 BIC=999.14
#
# forecasts
p2 <- forecast(fit2, 10, main = "ARIMA model")
p2$mean
# [1] 49636.48 50372.96 51109.45 51845.93 52582.41 53318.89 54055.37 54791.86
# [9] 55528.34 56264.82
Fabuła prognoz:
par(mfrow = c(2, 1), mar = c(2.5,2.2,2,2))
plot((cbind(x, p1$pred)), plot.type = "single", type = "n",
ylim = range(c(x, p1$pred + 1.96 * p1$se)), main = "Local trend model")
grid()
lines(x)
lines(p1$pred, col = "blue")
lines(p1$pred + 1.96 * p1$se, col = "red", lty = 2)
lines(p1$pred - 1.96 * p1$se, col = "red", lty = 2)
legend("topleft", legend = c("forecasts", "95% confidence interval"),
lty = c(1,2), col = c("blue", "red"), bty = "n")
plot((cbind(x, p2$mean)), plot.type = "single", type = "n",
ylim = range(c(x, p2$upper)), main = "ARIMA (0,1,0) with drift")
grid()
lines(x)
lines(p2$mean, col = "blue")
lines(ts(p2$lower[,2], start = end(x)[1] + 1), col = "red", lty = 2)
lines(ts(p2$upper[,2], start = end(x)[1] + 1), col = "red", lty = 2)
Prognozy są podobne w obu przypadkach i wyglądają rozsądnie. Zauważ, że prognozy są zgodne ze stosunkowo deterministycznym wzorcem podobnym do trendu liniowego, ale nie modelowaliśmy wyraźnie trendu liniowego. Powód jest następujący: i) w lokalnym modelu trendu wariancja składnika nachylenia jest szacowana na zero. To zmienia składową trendu w dryf, który daje efekt trendu liniowego. ii) ARIMA (0,1,1), model z dryfem jest wybierany w modelu dla szeregu różnicowego. Wpływ terminu stałego na szereg różnicowy jest trendem liniowym. Jest to omówione w tym poście .
Możesz sprawdzić, czy jeśli wybierzesz model lokalny lub ARIMA (0,1,0) bez dryfu, prognozy będą prostą linią poziomą, a zatem nie będą podobne do obserwowanej dynamiki danych. Cóż, jest to część układanki testów na pierwiastki jednostkowe i komponentów deterministycznych.
Edycja 1 (kontrola pozostałości):
autokorelacja i częściowy ACF nie sugerują struktury w pozostałościach.
resid1 <- residuals(fit1)
resid2 <- residuals(fit2)
par(mfrow = c(2, 2))
acf(resid1, lag.max = 20, main = "ACF residuals. Local trend model")
pacf(resid1, lag.max = 20, main = "PACF residuals. Local trend model")
acf(resid2, lag.max = 20, main = "ACF residuals. ARIMA(0,1,0) with drift")
pacf(resid2, lag.max = 20, main = "PACF residuals. ARIMA(0,1,0) with drift")
Jak sugeruje IrishStat, wskazane jest również sprawdzenie wartości odstających. Dwie dodatkowe wartości odstające są wykrywane za pomocą pakietu tsoutliers
.
require(tsoutliers)
resol <- tsoutliers(x, types = c("AO", "LS", "TC"),
remove.method = "bottom-up",
args.tsmethod = list(ic="bic", allowdrift=TRUE))
resol
#ARIMA(0,1,0) with drift
#Coefficients:
# drift AO2 AO51
# 736.4821 -3819.000 -4500.000
#s.e. 220.6171 1167.396 1167.397
#sigma^2 estimated as 2725622: log likelihood=-485.05
#AIC=978.09 AICc=978.88 BIC=986.2
#Outliers:
# type ind time coefhat tstat
#1 AO 2 2 -3819 -3.271
#2 AO 51 51 -4500 -3.855
Patrząc na ACF, możemy powiedzieć, że na poziomie 5% istotności, reszty są również losowe w tym modelu.
par(mfrow = c(2, 1))
acf(residuals(resol$fit), lag.max = 20, main = "ACF residuals. ARIMA with additive outliers")
pacf(residuals(resol$fit), lag.max = 20, main = "PACF residuals. ARIMA with additive outliers")
W tym przypadku obecność potencjalnych wartości odstających nie wydaje się zakłócać działania modeli. Potwierdza to test Jarque-Bera na normalność; zero normalności w resztach z początkowych modeli ( fit1
, fit2
) nie jest odrzucane na poziomie istotności 5%.
jarque.bera.test(resid1)[[1]]
# X-squared = 0.3221, df = 2, p-value = 0.8513
jarque.bera.test(resid2)[[1]]
#X-squared = 0.426, df = 2, p-value = 0.8082
Edycja 2 (wykres reszt i ich wartości)
Tak wyglądają resztki:
A oto ich wartości w formacie csv:
0;6.9205
-0.9571;-2942.4821
2.6108;4695.5179
-0.5453;-2065.4821
-0.2026;-771.4821
0.1242;-208.4821
0.1909;-120.4821
-0.0179;-535.4821
0.1449;-153.4821
0.484;526.5179
1.0748;1722.5179
0.3818;299.5179
-1.061;-2434.4821
0.0996;156.5179
0.4805;663.5179
0.8969;1463.5179
0.4111;461.5179
-1.0595;-2361.4821
0.0098;65.5179
0.5605;920.5179
0.8835;1481.5179
0.7669;1232.5179
1.4024;2593.5179
0.3785;473.5179
-1.1032;-2233.4821
-0.3813;-492.4821
2.2745;4642.5179
0.2935;154.5179
-0.1138;-165.4821
-0.8035;-1455.4821
-1.2982;-2321.4821
-1.9463;-3565.4821
-0.1648;62.5179
-0.1022;-253.4821
0.9755;1882.5179
-0.5662;-1418.4821
-0.0176;28.5179
0.5;928.5179
0.6831;1189.5179
-1.8889;-3964.4821
0.3896;1136.5179
-1.3113;-2799.4821
-0.9934;-1800.4821
-0.4085;-748.4821
1.2902;2482.5179
-0.0996;-657.4821
0.5539;981.5179
2.0007;3725.5179
1.0227;1490.5179
0.27;263.5179
-2.336;-4736.4821
1.8994;4263.5179
0.1301;-236.4821
-0.0892;-236.4821
-0.1148;-236.4821
-1.1207;-2236.4821
0.4801;1163.5179