library(readxl)
gdp <- read_excel("GDP per capita from 1800.xlsx")
View(gdp)
# загружаем файл GDP per capita from 1800
GDP<-log(gdp$Italy)
plot(GDP, type = "l")

acf(GDP, type='partial', na.action = na.omit)

Pacf(GDP)

ACF и PACF Значения графика ACF уходят в ноль после 1-ого лага,следовательно, параметр q процесса MA(q) равен 1. Также по графику автокорреляционной функции ACF заметно, что сезонность отсутствует.

Значения графика PACF уходят в ноль после первого лага, значит, параметр AR(p) процесса равен 1. Также по графику частной автокорреляционной функции для ВВП Италии, можно сказать, что значим только первый лаг, все остальные незначимы. Из найденных значений p и q получается, что модель ARIMA примет вид: ARIMA (1, d, 1).

Оценим параметры модели AR(1).

m<- Arima(GDP, c(1,0,0), include.constant =TRUE, method = c("CSS-ML"))  
coeftest(m)
## 
## z test of coefficients:
## 
##             Estimate Std. Error   z value  Pr(>|z|)    
## ar1       0.99930918 0.00093537 1068.3559 < 2.2e-16 ***
## intercept 9.10681470 1.20318526    7.5689 3.763e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
acf(residuals(m), na.action = na.omit)

Box.test(residuals(m), lag = 5, type = c("Ljung-Box"), fitdf = 1)
## 
##  Box-Ljung test
## 
## data:  residuals(m)
## X-squared = 33.531, df = 4, p-value = 9.3e-07

Остатки сильно коррелировны. Полученные оценки не являются несмещенными и состоятельными. P-value ниже 0.05, следовательно, нулевая гипотеза теста Льюинга-Бокса о том, что во временном ряду есть белый шум, не отвергается на уровне значимости 95%.

Тест Дики-Фуллера Для завершения построения модели ARIMA(p, d, q) необходимо определить порядок интеграции d. Для расчета этого значения проверим ряды на стационарность с помощью теста Дики-Фуллера. В случае обнаружения нестационарности будем рассматривать разности показателя. Количество взятых разностей для получения стационарного временного ряда и будет значением d - порядка интеграции.

#Спецификация модели без тренда(вокруг нуля)
ur.df(GDP, type="drift", lags=0)
## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 1.6041 8.8709
df<-ur.df(GDP, type="drift", lags=0)
summary(df)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.227690 -0.017150  0.000095  0.018405  0.287082 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.034792   0.029257  -1.189    0.236
## z.lag.1      0.005367   0.003346   1.604    0.110
## 
## Residual standard error: 0.04474 on 213 degrees of freedom
## Multiple R-squared:  0.01194,    Adjusted R-squared:  0.007298 
## F-statistic: 2.573 on 1 and 213 DF,  p-value: 0.1102
## 
## 
## Value of test-statistic is: 1.6041 8.8709 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
#Спецификация модели с трендом
ur.df(GDP, type="trend", lags=0)
## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -1.1997 7.6295 3.7237
df<-ur.df(GDP, type="trend", lags=0)
summary(df)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.240460 -0.014509  0.000978  0.020015  0.269262 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.0616939  0.0526227   1.172   0.2424  
## z.lag.1     -0.0086023  0.0071705  -1.200   0.2316  
## tt           0.0002315  0.0001054   2.197   0.0291 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04435 on 212 degrees of freedom
## Multiple R-squared:  0.03394,    Adjusted R-squared:  0.02482 
## F-statistic: 3.724 on 2 and 212 DF,  p-value: 0.02574
## 
## 
## Value of test-statistic is: -1.1997 7.6295 3.7237 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

Видно, что гипотеза о нестационарности не может быть отвегнута(-1.199 >-3.43).

# расширенный тест Дики-Фуллера, KPSS и PP тесты.
ur.df(GDP, type="drift", lags=2)
## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 1.0313 4.2975
df<-ur.df(GDP, type="drift", lags=2)
summary(df)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.161165 -0.017400  0.000049  0.017372  0.307683 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.020181   0.027678  -0.729   0.4667    
## z.lag.1      0.003283   0.003183   1.031   0.3036    
## z.diff.lag1  0.423035   0.068727   6.155 3.77e-09 ***
## z.diff.lag2 -0.130836   0.069014  -1.896   0.0594 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0415 on 209 degrees of freedom
## Multiple R-squared:  0.1641, Adjusted R-squared:  0.1521 
## F-statistic: 13.67 on 3 and 209 DF,  p-value: 3.54e-08
## 
## 
## Value of test-statistic is: 1.0313 4.2975 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
kpss.test (GDP, null="Level")
## Warning in kpss.test(GDP, null = "Level"): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  GDP
## KPSS Level = 3.6622, Truncation lag parameter = 4, p-value = 0.01
kpss.test (GDP, null="Trend")
## Warning in kpss.test(GDP, null = "Trend"): p-value smaller than printed p-value
## 
##  KPSS Test for Trend Stationarity
## 
## data:  GDP
## KPSS Trend = 0.98824, Truncation lag parameter = 4, p-value = 0.01
pp.test(GDP)
## 
##  Phillips-Perron Unit Root Test
## 
## data:  GDP
## Dickey-Fuller Z(alpha) = -2.5091, Truncation lag parameter = 4, p-value
## = 0.9536
## alternative hypothesis: stationary
Acf(GDP)

Pacf(GDP)

Таким образом, все тесты не противоречат друг другу и указывают на нестационарность ряда (которую видно визуально).

# переходим к темпам прироста ВВП на душу
dGDP<-diff(GDP, differences=2)
plot(dGDP, type = "l")

ur.df(dGDP, type="drift", lags=0)
## 
## ############################################################### 
## # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -18.2799 167.0788
df1<-ur.df(dGDP, type="drift", lags=0)
summary(df1)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12805 -0.01780 -0.00127  0.01733  0.43037 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0001059  0.0033492  -0.032    0.975    
## z.lag.1     -1.2242845  0.0669742 -18.280   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04888 on 211 degrees of freedom
## Multiple R-squared:  0.613,  Adjusted R-squared:  0.6111 
## F-statistic: 334.2 on 1 and 211 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -18.2799 167.0788 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
kpss.test (dGDP, null="Level")
## Warning in kpss.test(dGDP, null = "Level"): p-value greater than printed
## p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  dGDP
## KPSS Level = 0.012893, Truncation lag parameter = 4, p-value = 0.1
pp.test(dGDP)
## Warning in pp.test(dGDP): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  dGDP
## Dickey-Fuller Z(alpha) = -200.54, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
Acf(dGDP)

Pacf(dGDP)

Видно, что значимы 1-6,8 и 10 лаги. Шестнадцатый лаг достаточно отдаленный, имеет смысл построить модель AR порядка 6,8 и 10.

#Оценим параметры модели AR(6)
m6 <- Arima(dGDP, c(6,0,0), include.constant =TRUE, method = c("CSS-ML"))  
coeftest(m6) 
## 
## z test of coefficients:
## 
##              Estimate  Std. Error z value  Pr(>|z|)    
## ar1       -4.5156e-01  6.7657e-02 -6.6743 2.484e-11 ***
## ar2       -4.9842e-01  7.2837e-02 -6.8429 7.759e-12 ***
## ar3       -4.3606e-01  7.7817e-02 -5.6037 2.098e-08 ***
## ar4       -2.8574e-01  7.7848e-02 -3.6705 0.0002421 ***
## ar5       -2.1117e-01  7.2485e-02 -2.9133 0.0035765 ** 
## ar6       -1.4940e-01  6.7643e-02 -2.2087 0.0271973 *  
## intercept -9.1655e-05  9.8915e-04 -0.0927 0.9261732    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Темп прироста ВВП предыдущих 6 лет отрицательно влияют на темп прироста ВВП текущего года на уровне значимости 0.05.

#Оценим параметры модели AR(8)
m8 <- Arima(dGDP, c(8,0,0), include.constant =TRUE, method = c("CSS-ML"))  
coeftest(m8)
## 
## z test of coefficients:
## 
##              Estimate  Std. Error z value  Pr(>|z|)    
## ar1       -0.49290109  0.06585383 -7.4848 7.167e-14 ***
## ar2       -0.57300264  0.07236779 -7.9179 2.415e-15 ***
## ar3       -0.53563333  0.07900681 -6.7796 1.205e-11 ***
## ar4       -0.41924145  0.08307506 -5.0465 4.499e-07 ***
## ar5       -0.38199564  0.08254565 -4.6277 3.698e-06 ***
## ar6       -0.33202978  0.07872055 -4.2178 2.467e-05 ***
## ar7       -0.21857076  0.07194683 -3.0379  0.002382 ** 
## ar8       -0.26616372  0.06577061 -4.0468 5.191e-05 ***
## intercept -0.00010817  0.00068606 -0.1577  0.874723    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Темпы прироста за последние 8 лет сдерживают темп прироста ВВП в текущем периоде на уровне значимости 0.01.

#Оценим параметры модели AR(10)
m10 <- Arima(dGDP, c(10,0,0), include.constant =TRUE, method = c("CSS-ML"))  
coeftest(m10)
## 
## z test of coefficients:
## 
##              Estimate  Std. Error z value  Pr(>|z|)    
## ar1       -4.9519e-01  6.7323e-02 -7.3554 1.904e-13 ***
## ar2       -6.2039e-01  7.5126e-02 -8.2580 < 2.2e-16 ***
## ar3       -5.7644e-01  8.2582e-02 -6.9802 2.947e-12 ***
## ar4       -4.8003e-01  8.9335e-02 -5.3734 7.726e-08 ***
## ar5       -4.5318e-01  9.0656e-02 -4.9989 5.764e-07 ***
## ar6       -4.0790e-01  9.0203e-02 -4.5221 6.123e-06 ***
## ar7       -3.1232e-01  8.8339e-02 -3.5354 0.0004071 ***
## ar8       -3.6677e-01  8.2256e-02 -4.4588 8.241e-06 ***
## ar9       -8.8725e-02  7.5465e-02 -1.1757 0.2397099    
## ar10      -1.7095e-01  6.7668e-02 -2.5263 0.0115271 *  
## intercept -9.2612e-05  5.7630e-04 -0.1607 0.8723297    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Темпы прироста за последние 10 лет сдерживают темп прироста ВВП в текущем периоде.

#Проверим автокорреляции остатков
acf(residuals(m6), na.action = na.omit)

Box.test(residuals(m6), lag = 8, type = c("Ljung-Box"), fitdf = 6)
## 
##  Box-Ljung test
## 
## data:  residuals(m6)
## X-squared = 16.488, df = 2, p-value = 0.0002628
acf(residuals(m8), na.action = na.omit)

Box.test(residuals(m8), lag = 0, type = c("Ljung-Box"), fitdf = 8)
## 
##  Box-Ljung test
## 
## data:  residuals(m8)
## X-squared = NA, df = -8, p-value = NA
acf(residuals(m10), na.action = na.omit)

Box.test(residuals(m10), lag = 0, type = c("Ljung-Box"), fitdf = 10)
## 
##  Box-Ljung test
## 
## data:  residuals(m10)
## X-squared = NA, df = -10, p-value = NA

ACF и тесты Льюнга- Бокса показывают, что гипотезу об отсутствии автокорреляция остатков для 8 лагов можно отвергнуть на уровне значимости 0.1%, следовательно, остатки не распределены независимо и коррелируют.

# создаем лаги для оценивания "коротких" моделей
dGDP1 <- c(0,dGDP[1:length(dGDP)-1])
dGDP2 <- c(0,0,dGDP[2:length(dGDP)-2])
dGDP3 <- c(0,0,0,dGDP[3:length(dGDP)-3])
dGDP4 <- c(0,0,0,0,dGDP[4:length(dGDP)-4])
dGDP5<-  c(0,0,0,0,0,dGDP[5:length(dGDP)-5])
dGDP6<-  c(0,0,0,0,0,0,dGDP[6:length(dGDP)-6])
dGDP8<-  c(0,0,0,0,0,0,0,0,dGDP[8:length(dGDP)-8])
m8short=lm(dGDP ~ dGDP1 + dGDP2 + dGDP3 + dGDP4 + dGDP5 + dGDP6 + dGDP8)
summary(m8short)
## 
## Call:
## lm(formula = dGDP ~ dGDP1 + dGDP2 + dGDP3 + dGDP4 + dGDP5 + dGDP6 + 
##     dGDP8)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16129 -0.01750 -0.00137  0.01648  0.34130 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.979e-05  2.965e-03   0.020 0.983930    
## dGDP1       -4.559e-01  6.756e-02  -6.748 1.49e-10 ***
## dGDP2       -5.165e-01  7.302e-02  -7.073 2.32e-11 ***
## dGDP3       -4.620e-01  7.835e-02  -5.896 1.50e-08 ***
## dGDP4       -3.121e-01  7.836e-02  -3.983 9.45e-05 ***
## dGDP5       -2.573e-01  7.437e-02  -3.460 0.000658 ***
## dGDP6       -2.073e-01  7.058e-02  -2.936 0.003697 ** 
## dGDP8       -1.827e-01  6.270e-02  -2.913 0.003969 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04337 on 206 degrees of freedom
## Multiple R-squared:  0.2729, Adjusted R-squared:  0.2482 
## F-statistic: 11.04 on 7 and 206 DF,  p-value: 7.814e-12
acf(residuals(m8short), na.action = na.omit)

Box.test(residuals(m8short), lag = 7, type = c("Ljung-Box"), fitdf = 7)
## 
##  Box-Ljung test
## 
## data:  residuals(m8short)
## X-squared = 11.621, df = 0, p-value < 2.2e-16

Гипотезу о равенстве автокорреляции нулю для короткой модели отвергнуть нельзя.

AIC(m6,m8,m10,m8short)
##         df       AIC
## m6       8 -718.9306
## m8      10 -732.6910
## m10     12 -734.9668
## m8short  9 -725.8889
BIC(m6,m8,m10,m8short)
##         df       BIC
## m6       8 -692.0028
## m8      10 -699.0312
## m10     12 -694.5751
## m8short  9 -695.5951