Экзаменационный проект
Хасенова Снежана, Малиновская Ирина
Загружаем библиотеки и данные:
library('tseries')
library(readxl)
library("vars")
library('TSA')
library(lmtest)
library("portes")
library('urca')
library('forecast')
library("portes")
library(dplyr)
library("strucchange")
library("vars")
data= read_excel("data_LAI.xlsx") В нашем датасете всё объединено, разделим по странам в разные датафреймы:
argentina= data %>% filter(data$countrycode== 'ARG')
bolivia= data %>% filter(countrycode== 'BOL')
brazil= data %>% filter(countrycode== 'BRA')
colombia= data %>% filter(countrycode== 'COL')
cuba= data %>% filter(countrycode== 'CUB')Построим графики ВВП на душу населения:
plot(argentina$year, argentina$gdppc, col = 'red', type = 'l', xlab = 'Год', ylab = 'ВВП, $', main = 'ВВП на душу населения, Аргентина')
plot(bolivia$year, bolivia$gdppc, col = 'orange', type = 'l', xlab = 'Год', ylab = 'ВВП, $', main = 'ВВП на душу населения, Боливия')
plot(brazil$year, brazil$gdppc, col = 'green', type = 'l', xlab = 'Год', ylab = 'ВВП, $', main = 'ВВП на душу населения, Бразилия')
plot(colombia$year, colombia$gdppc, col = 'blue', type = 'l', xlab = 'Год', ylab = 'ВВП, $', main = 'ВВП на душу населения, Колумбия')
plot(cuba$year, cuba$gdppc, col = 'purple', type = 'l', xlab = 'Год', ylab = 'ВВП, $', main = 'ВВП на душу населения, Куба')
### Построим графики динамики роста населения:
plot(argentina$year, argentina$pop, col='red', type = 'l', xlab='Год', ylab='Население, тыс.',
main = 'Динамика роста населения, Аргентина')
plot(bolivia$year, bolivia$pop, col='orange', type = 'l', xlab='Год', ylab='Население, тыс.',
main = 'Динамика роста населения, Боливия')
plot(brazil$year, brazil$pop, col='green', type = 'l', xlab='Год', ylab='Население, тыс.',
main = 'Динамика роста населения, Бразилия')
plot(colombia$year, colombia$pop, col='blue', type = 'l', xlab='Год', ylab='Население, тыс.',
main = 'Динамика роста населения, Колумбия')
plot(cuba$year, cuba$pop, col='purple', type = 'l',xlab='Год', ylab='Население, тыс.',
main = 'Динамика роста населения, Куба')
Прологарифмируем динамику роста ВВП и рассчитаем первые разности:
argentina= argentina %>% mutate(log_dyn= log(dynamics_gdp))
d1arg= diff(argentina$log_dyn, differences = 1)
bolivia= bolivia %>% mutate(log_dyn= log(dynamics_gdp))
d1bol= diff(bolivia$log_dyn, differences = 1)
brazil= brazil %>% mutate(log_dyn= log(dynamics_gdp))
d1bra= diff(brazil$log_dyn, differences = 1)
colombia= colombia %>% mutate(log_dyn= log(dynamics_gdp))
d1col= diff(colombia$log_dyn, differences = 1)
cuba= cuba %>% mutate(log_dyn= log(dynamics_gdp))
d1cub= diff(cuba$log_dyn, differences = 1)Проверим стационарность:
Аргентина
pacf(argentina$log_dyn)
st_arg= ur.df(argentina$log_dyn, type='drift', lags=1, selectlags = 'Fixed')
summary(st_arg)##
## ###############################################
## # 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.12469 -0.04016 0.01116 0.03479 0.09164
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.176523 0.115677 1.526 0.130
## z.lag.1 -0.007853 0.005998 -1.309 0.194
## z.diff.lag 0.125162 0.101094 1.238 0.219
##
## Residual standard error: 0.04952 on 96 degrees of freedom
## Multiple R-squared: 0.03493, Adjusted R-squared: 0.01482
## F-statistic: 1.737 on 2 and 96 DF, p-value: 0.1815
##
##
## Value of test-statistic is: -1.3092 10.2906
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Боливия
pacf(bolivia$log_dyn)
st_bol= ur.df(bolivia$log_dyn, type='drift', lags=1, selectlags = 'Fixed')
summary(st_bol)##
## ###############################################
## # 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.130933 -0.012211 0.006898 0.017859 0.094601
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.054089 0.079121 -0.684 0.49586
## z.lag.1 0.004556 0.004827 0.944 0.34760
## z.diff.lag 0.289574 0.098602 2.937 0.00415 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03804 on 96 degrees of freedom
## Multiple R-squared: 0.1034, Adjusted R-squared: 0.08467
## F-statistic: 5.533 on 2 and 96 DF, p-value: 0.00532
##
##
## Value of test-statistic is: 0.9439 9.7206
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Бразилия
pacf(brazil$log_dyn)
st_bra= ur.df(brazil$log_dyn, type='drift', lags=1, selectlags = 'Fixed')
summary(st_bra)##
## ###############################################
## # 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.162851 -0.020930 0.004812 0.026798 0.102004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.091593 0.055059 1.664 0.09947 .
## z.lag.1 -0.002949 0.002751 -1.072 0.28640
## z.diff.lag 0.265111 0.096811 2.738 0.00736 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04067 on 96 degrees of freedom
## Multiple R-squared: 0.09062, Adjusted R-squared: 0.07168
## F-statistic: 4.783 on 2 and 96 DF, p-value: 0.01046
##
##
## Value of test-statistic is: -1.072 15.1377
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Колумбия
pacf(colombia$log_dyn)
st_col= ur.df(colombia$log_dyn, type='drift', lags=1, selectlags = 'Fixed')
summary(st_col)##
## ###############################################
## # 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.071697 -0.013686 0.004168 0.012029 0.075156
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.092456 0.036990 2.499 0.0141 *
## z.lag.1 -0.003324 0.001951 -1.704 0.0917 .
## z.diff.lag 0.251546 0.097406 2.582 0.0113 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02309 on 96 degrees of freedom
## Multiple R-squared: 0.1122, Adjusted R-squared: 0.09375
## F-statistic: 6.069 on 2 and 96 DF, p-value: 0.003296
##
##
## Value of test-statistic is: -1.7038 21.8179
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Куба
pacf(cuba$log_dyn)
st_cub= ur.df(cuba$log_dyn, type='drift', lags=1, selectlags = 'Fixed')
summary(st_cub)##
## ###############################################
## # 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.293576 -0.025075 0.004754 0.035272 0.112526
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.067199 0.136553 0.492 0.623765
## z.lag.1 -0.002954 0.008043 -0.367 0.714253
## z.diff.lag 0.343171 0.095922 3.578 0.000546 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06287 on 96 degrees of freedom
## Multiple R-squared: 0.1178, Adjusted R-squared: 0.09942
## F-statistic: 6.41 on 2 and 96 DF, p-value: 0.002439
##
##
## Value of test-statistic is: -0.3672 3.223
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81Проверим стационарность первых разностей:
Аргентина
pacf(d1arg)
d1_arg= ur.df(d1arg, type='drift', lags=8, selectlags = 'Fixed')
summary(d1_arg)##
## ###############################################
## # 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.106229 -0.034031 0.003702 0.032841 0.078131
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.03572 0.01084 3.296 0.001457 **
## z.lag.1 -1.28539 0.32682 -3.933 0.000176 ***
## z.diff.lag1 0.52886 0.29271 1.807 0.074507 .
## z.diff.lag2 0.35018 0.27963 1.252 0.214063
## z.diff.lag3 0.35893 0.24691 1.454 0.149885
## z.diff.lag4 0.28117 0.22123 1.271 0.207403
## z.diff.lag5 0.24929 0.19446 1.282 0.203515
## z.diff.lag6 0.07968 0.16896 0.472 0.638498
## z.diff.lag7 0.35007 0.13287 2.635 0.010090 *
## z.diff.lag8 -0.05731 0.11170 -0.513 0.609305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04672 on 81 degrees of freedom
## Multiple R-squared: 0.5615, Adjusted R-squared: 0.5128
## F-statistic: 11.53 on 9 and 81 DF, p-value: 1.962e-11
##
##
## Value of test-statistic is: -3.9331 7.842
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86Боливия
pacf(d1bol)
d1_bol= ur.df(d1bol, type='drift', lags=3, selectlags = 'Fixed')
summary(d1_bol)##
## ###############################################
## # 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.112909 -0.015606 0.005949 0.016523 0.080524
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.014743 0.005267 2.799 0.006252 **
## z.lag.1 -0.502325 0.134027 -3.748 0.000313 ***
## z.diff.lag1 -0.081590 0.136939 -0.596 0.552780
## z.diff.lag2 -0.243812 0.113341 -2.151 0.034115 *
## z.diff.lag3 0.187675 0.098545 1.904 0.060010 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03425 on 91 degrees of freedom
## Multiple R-squared: 0.4602, Adjusted R-squared: 0.4364
## F-statistic: 19.39 on 4 and 91 DF, p-value: 1.442e-11
##
##
## Value of test-statistic is: -3.7479 7.0235
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86Бразилия
pacf(d1bra)
d1_bra= ur.df(d1bra, type='drift', lags=1, selectlags = 'Fixed')
summary(d1_bra)##
## ###############################################
## # 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.161480 -0.021098 0.002037 0.028047 0.099340
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.026926 0.006966 3.865 0.000203 ***
## z.lag.1 -0.602445 0.120291 -5.008 2.53e-06 ***
## z.diff.lag -0.182482 0.099400 -1.836 0.069508 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04032 on 95 degrees of freedom
## Multiple R-squared: 0.3902, Adjusted R-squared: 0.3774
## F-statistic: 30.4 on 2 and 95 DF, p-value: 6.239e-11
##
##
## Value of test-statistic is: -5.0082 12.5978
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86Колумбия
pacf(d1col)
d1_col= ur.df(d1col, type='drift', lags=1, selectlags = 'Fixed')
summary(d1_col)##
## ###############################################
## # 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.074975 -0.014949 0.000419 0.014117 0.081128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.029537 0.005734 5.151 1.40e-06 ***
## z.lag.1 -0.706651 0.121773 -5.803 8.55e-08 ***
## z.diff.lag -0.025372 0.101514 -0.250 0.803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02348 on 95 degrees of freedom
## Multiple R-squared: 0.3654, Adjusted R-squared: 0.3521
## F-statistic: 27.35 on 2 and 95 DF, p-value: 4.153e-10
##
##
## Value of test-statistic is: -5.803 16.8887
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86Куба
pacf(d1cub)
d1_cub= ur.df(d1cub, type='drift', lags=1, selectlags = 'Fixed')
summary(d1_cub)##
## ###############################################
## # 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.291244 -0.025167 0.003958 0.037338 0.111148
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.017137 0.007108 2.411 0.0178 *
## z.lag.1 -0.656339 0.117671 -5.578 2.29e-07 ***
## z.diff.lag -0.002961 0.102455 -0.029 0.9770
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06323 on 95 degrees of freedom
## Multiple R-squared: 0.3291, Adjusted R-squared: 0.315
## F-statistic: 23.3 on 2 and 95 DF, p-value: 5.826e-09
##
##
## Value of test-statistic is: -5.5777 15.5564
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86ARMA- модели для первых разностей лог(ввп*популярность):
Аргентина
eacf(d1arg)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o x o o o o o o
## 1 x o o o o o o x o o o o o o
## 2 x o o o o o o x o o o o o o
## 3 x o x o o o o x o o o o o o
## 4 x x o o o o o x o o o o o o
## 5 x x o x o o o x o o o o o o
## 6 x o o x o o o o o o o o o o
## 7 x o x x o o o o o o o o o o#ARMA (0,0)
arma_arg00= Arima(d1arg, c(0,1,0), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_arg00) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## drift -0.00061662 0.00654973 -0.0941 0.925summary(arma_arg00) ## Series: d1arg
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## -0.0006
## s.e. 0.0065
##
## sigma^2 = 0.004289: log likelihood = 129.88
## AIC=-255.77 AICc=-255.64 BIC=-250.58
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.652884e-07 0.06483469 0.05114007 148.6809 282.0843 0.9889339
## ACF1
## Training set -0.3624866Acf(residuals(arma_arg00)) 
Box.test(residuals(arma_arg00), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 0.0007433##
## Box-Ljung test
##
## data: residuals(arma_arg00)
## X-squared = 19.123, df = 4, p-value = 0.0007433#ARMA (0,1)
arma_arg0= Arima(d1arg, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_arg0) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.99999948 0.02777698 -36.0010 <2e-16 ***
## drift -0.00018792 0.00019836 -0.9474 0.3434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_arg0) ## Series: d1arg
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## -1.0000 -2e-04
## s.e. 0.0278 2e-04
##
## sigma^2 = 0.002485: log likelihood = 155.12
## AIC=-304.23 AICc=-303.98 BIC=-296.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001337162 0.04909204 0.04020988 64.6865 169.8341 0.7775688
## ACF1
## Training set 0.1251532Acf(residuals(arma_arg0)) 
Box.test(residuals(arma_arg0), lag = 5, type = c("Ljung-Box"), fitdf = 1)#p-value = 0.4084##
## Box-Ljung test
##
## data: residuals(arma_arg0)
## X-squared = 3.9823, df = 4, p-value = 0.4084Acf(residuals(arma_arg0)^2)
Box.test(residuals(arma_arg0)^2, lag = 5, type = c("Ljung-Box"), fitdf =2) #p-value = 0.0936##
## Box-Ljung test
##
## data: residuals(arma_arg0)^2
## X-squared = 4.7562, df = 3, p-value = 0.1905#SUP-F
supf_arg= Fstats(d1arg ~ 1, from = 0.1, to = NULL)
plot(supf_arg , alpha = 0.01) + lines(breakpoints(supf_arg ))
## integer(0)sup= breakpoints(supf_arg )
sup$breakpoints #11## [1] 11sctest(supf_arg , type = "supF") #p-value = 0.5515##
## supF test
##
## data: supf_arg
## sup.F = 3.3623, p-value = 0.5515Прогноз:
forecast_arg0= forecast(arma_arg0, h=30)
plot(forecast_arg0, col='red')
#ARMA модель на периоде (1,11)
arma_arg_11 = d1arg[seq(1,11)]
sup_11 = Fstats(arma_arg_11 ~ 1, from = 0.1, to = NULL)## Warning in Fstats(arma_arg_11 ~ 1, from = 0.1, to = NULL): 'from' changed (was
## too small)## Warning in Fstats(arma_arg_11 ~ 1, from = 0.1, to = NULL): 'to' changed (was too
## large)plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #6## [1] 6sctest(sup_11 , type = "supF") #p-value = 0.8796##
## supF test
##
## data: sup_11
## sup.F = 1.4883, p-value = 0.8796#ARMA модель на периоде (11,60)
arma_arg_60 = d1arg[seq(11,60)]
sup_11 = Fstats(arma_arg_60 ~ 1, from = 0.1, to = NULL)
plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #5## [1] 5sctest(sup_11 , type = "supF") #p-value = 0.296##
## supF test
##
## data: sup_11
## sup.F = 4.9742, p-value = 0.296#ARMA модель на периоде (60,100)
arma_arg_100 = d1arg[seq(60,100)]
sup_11 = Fstats(arma_arg_100 ~ 1, from = 0.1, to = NULL)
plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #13## [1] 13sctest(sup_11 , type = "supF") #p-value = 0.4731##
## supF test
##
## data: sup_11
## sup.F = 3.8, p-value = 0.4731Боливия
eacf(d1bol)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x o o o o o o o o o o o
## 1 o o x o o o o o o o o o o o
## 2 o x x o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x o o o o o o o o o o o o o
## 5 x o o o o o o o o o o o o o
## 6 x o o o o x o o x o o o o o
## 7 x o o o x x o o o o o o o o#ARMA (0,1)
arma_bol0= Arima(d1bol, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_bol0) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.72121208 0.09832614 -7.3349 2.219e-13 ***
## drift 0.00029295 0.00113760 0.2575 0.7968
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bol0) ## Series: d1bol
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## -0.7212 0.0003
## s.e. 0.0983 0.0011
##
## sigma^2 = 0.001582: log likelihood = 179.4
## AIC=-352.8 AICc=-352.55 BIC=-345.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.000170826 0.0391713 0.02731505 321.2612 357.422 0.8916132
## ACF1
## Training set 0.08634149Acf(residuals(arma_bol0)) 
Box.test(residuals(arma_bol0), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 0.001262##
## Box-Ljung test
##
## data: residuals(arma_bol0)
## X-squared = 17.951, df = 4, p-value = 0.001262#ARMA (1,0)
arma_bol1= Arima(d1bol, c(1,1,0), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_bol1) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.33034362 0.09420528 -3.5066 0.0004538 ***
## drift 0.00030071 0.00332450 0.0905 0.9279264
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bol1) ## Series: d1bol
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## -0.3303 0.0003
## s.e. 0.0942 0.0033
##
## sigma^2 = 0.001965: log likelihood = 168.98
## AIC=-331.96 AICc=-331.71 BIC=-324.18
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -2.284134e-05 0.04365505 0.02793162 238.3325 286.6156 0.9117392
## ACF1
## Training set -0.1807575Acf(residuals(arma_bol1)) 
Box.test(residuals(arma_bol1), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 5.035e-09##
## Box-Ljung test
##
## data: residuals(arma_bol1)
## X-squared = 44.507, df = 4, p-value = 5.035e-09#ARMA (1,1)
arma_bol11= Arima(d1bol, c(1,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_bol11) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.29785153 0.09653826 3.0853 0.002033 **
## ma1 -0.99999726 0.02881228 -34.7073 < 2.2e-16 ***
## drift 0.00023687 0.00020791 1.1393 0.254566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bol11) ## Series: d1bol
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.2979 -1.0000 2e-04
## s.e. 0.0965 0.0288 2e-04
##
## sigma^2 = 0.001439: log likelihood = 182.98
## AIC=-357.96 AICc=-357.54 BIC=-347.58
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0003003473 0.03716281 0.02541496 248.2704 278.314 0.8295908
## ACF1
## Training set 0.002607564Acf(residuals(arma_bol11)) 
Box.test(residuals(arma_bol11), lag = 5, type = c("Ljung-Box"), fitdf = 2) #p-value = 0.0001222##
## Box-Ljung test
##
## data: residuals(arma_bol11)
## X-squared = 20.688, df = 3, p-value = 0.0001222#ARMA (2,0)
arma_bol2= Arima(d1bol, c(2,1,0), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_bol2) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.51043618 0.08443525 -6.0453 1.491e-09 ***
## ar2 -0.53201926 0.08354324 -6.3682 1.913e-10 ***
## drift 0.00042385 0.00183408 0.2311 0.8172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bol2) ## Series: d1bol
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## -0.5104 -0.5320 0.0004
## s.e. 0.0844 0.0835 0.0018
##
## sigma^2 = 0.001406: log likelihood = 185.72
## AIC=-363.45 AICc=-363.02 BIC=-353.06
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.000197172 0.03673942 0.02558803 330.7571 391.3236 0.83524
## ACF1
## Training set 0.02969799Acf(residuals(arma_bol2)) 
Box.test(residuals(arma_bol2), lag = 5, type = c("Ljung-Box"), fitdf = 2) #p-value = 0.1564##
## Box-Ljung test
##
## data: residuals(arma_bol2)
## X-squared = 5.219, df = 3, p-value = 0.1564Acf(residuals(arma_bol2)^2)
Box.test(residuals(arma_bol2)^2, lag = 5, type = c("Ljung-Box"), fitdf =2) #p-value = 1.338e-07##
## Box-Ljung test
##
## data: residuals(arma_bol2)^2
## X-squared = 34.807, df = 3, p-value = 1.338e-07Прогноз:
forecast_bol2= forecast(arma_bol2, h=30)
plot(forecast_bol2, col='red')
#SUP-F
supf_bol= Fstats(d1bol ~ 1, from = 0.1, to = NULL)
plot(supf_bol , alpha = 0.01) + lines(breakpoints(supf_bol ))
## integer(0)sup= breakpoints(supf_bol )
sup$breakpoints #41## [1] 41sctest(supf_bol , type = "supF") #p-value = 0.1892##
## supF test
##
## data: supf_bol
## sup.F = 6.0441, p-value = 0.1892#ARMA модель на периоде (1,41)
arma_bol_41 = d1bol[seq(1,41)]
sup_bol = Fstats(arma_bol_41 ~ 1, from = 0.1, to = NULL)
plot(sup_bol , alpha = 0.01) + lines(breakpoints(sup_bol ))
## integer(0)breakpoints(sup_bol )$breakpoints #31## [1] 31sctest(sup_bol , type = "supF") #p-value = 0.4838##
## supF test
##
## data: sup_bol
## sup.F = 3.7403, p-value = 0.4838#ARMA модель на периоде (41,60)
arma_bol_60 = d1bol[seq(41,60)]
sup_bol = Fstats(arma_bol_60 ~ 1, from = 0.1, to = NULL)
plot(sup_bol , alpha = 0.01) + lines(breakpoints(sup_bol ))
## integer(0)breakpoints(sup_bol )$breakpoints #3## [1] 3sctest(sup_bol , type = "supF") #p-value = 0.002647##
## supF test
##
## data: sup_bol
## sup.F = 15.25, p-value = 0.002647#ARMA модель на периоде (60,100)
arma_bol_100 = d1bol[seq(60,100)]
sup_bol = Fstats(arma_bol_100 ~ 1, from = 0.1, to = NULL)
plot(sup_bol , alpha = 0.01) + lines(breakpoints(sup_bol ))
## integer(0)breakpoints(sup_bol )$breakpoints #9## [1] 9sctest(sup_bol , type = "supF") #p-value = 2.284e-14##
## supF test
##
## data: sup_bol
## sup.F = 65.462, p-value = 2.284e-14Бразилия
eacf(d1bra)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o o#ARMA (1,1)
arma_bra11= Arima(d1bra, c(1,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_bra11) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.27993477 0.09854763 2.8406 0.004503 **
## ma1 -0.99999933 0.09529360 -10.4939 < 2.2e-16 ***
## drift -0.00027070 0.00021788 -1.2424 0.214079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bra11) ## Series: d1bra
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.2799 -1.0000 -3e-04
## s.e. 0.0985 0.0953 2e-04
##
## sigma^2 = 0.001698: log likelihood = 174.76
## AIC=-341.52 AICc=-341.1 BIC=-331.14
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0002025987 0.04037175 0.03049277 -1524.666 1679.205 0.8030053
## ACF1
## Training set -0.04598922Acf(residuals(arma_bra11)) 
Box.test(residuals(arma_bra11), lag = 5, type = c("Ljung-Box"), fitdf = 2) #p-value = 0.1223##
## Box-Ljung test
##
## data: residuals(arma_bra11)
## X-squared = 5.7902, df = 3, p-value = 0.1223Acf(residuals(arma_bra11)^2)
Box.test(residuals(arma_bra11)^2, lag = 5, type = c("Ljung-Box"), fitdf =2) #p-value = 0.8865##
## Box-Ljung test
##
## data: residuals(arma_bra11)^2
## X-squared = 0.6431, df = 3, p-value = 0.8865Прогноз:
forecast_bra11= forecast(arma_bra11, h=30)
plot(forecast_bra11, col='red')
#SUP-F
supf_bra= Fstats(d1bra ~ 1, from = 0.1, to = NULL)
plot(supf_bra , alpha = 0.01) + lines(breakpoints(supf_bra ))
## integer(0)sup= breakpoints(supf_bra )
sup$breakpoints #62## [1] 62sctest(supf_bra , type = "supF") #p-value = 0.007449##
## supF test
##
## data: supf_bra
## sup.F = 13.1, p-value = 0.007449#ARMA модель на периоде (1,62)
arma_bra_62 = d1bra[seq(1,62)]
supf_bra = Fstats(arma_bra_62 ~ 1, from = 0.1, to = NULL)
plot(supf_bra , alpha = 0.01) + lines(breakpoints(supf_bra ))
## integer(0)breakpoints(supf_bra )$breakpoints #24## [1] 24sctest(supf_bra , type = "supF") #p-value = 0.01689##
## supF test
##
## data: supf_bra
## sup.F = 11.398, p-value = 0.01689#ARMA модель на периоде (62,100)
arma_bra_100 = d1bra[seq(62,100)]
supf_bra = Fstats(arma_bra_100 ~ 1, from = 0.1, to = NULL)
plot(supf_bra , alpha = 0.01) + lines(breakpoints(supf_bra ))
## integer(0)breakpoints(supf_bra )$breakpoints #36## [1] 36sctest(supf_bra , type = "supF") #p-value = 0.04773##
## supF test
##
## data: supf_bra
## sup.F = 9.3686, p-value = 0.04773Колумбия
eacf(d1col)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 x o o o o o o o o o o o o o
## 3 x o x o o o o o o o o o o o
## 4 o x x x o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x o o o o o o o o o o o o o
## 7 x o o o o o o o o o o o o o#ARMA (0,1)
arma_col0= Arima(d1col, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_col0) ## Warning in sqrt(diag(se)): NaNs produced##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.91686155 NaN NaN NaN
## drift -0.00026920 0.00011345 -2.3728 0.01765 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_col0) ## Series: d1col
## ARIMA(0,1,1) with drift
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ma1 drift
## -0.9169 -3e-04
## s.e. NaN 1e-04
##
## sigma^2 = 0.0005989: log likelihood = 226.93
## AIC=-447.86 AICc=-447.61 BIC=-440.07
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.00123263 0.02410197 0.01783244 -24.0679 96.9182 0.8572924
## ACF1
## Training set 0.2277591Acf(residuals(arma_col0)) 
Box.test(residuals(arma_col0), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 0.02433##
## Box-Ljung test
##
## data: residuals(arma_col0)
## X-squared = 11.207, df = 4, p-value = 0.02433#ARMA (1,0)
arma_col1= Arima(d1col, c(1,1,0), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_col1) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.3717111 0.0927607 -4.0072 6.144e-05 ***
## drift -0.0006076 0.0019707 -0.3083 0.7578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_col1) ## Series: d1col
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## -0.3717 -6e-04
## s.e. 0.0928 2e-03
##
## sigma^2 = 0.0007324: log likelihood = 217.81
## AIC=-429.62 AICc=-429.37 BIC=-421.84
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.537389e-05 0.02665369 0.01955593 -10.05598 72.55121 0.9401489
## ACF1
## Training set -0.1177631Acf(residuals(arma_col1)) 
Box.test(residuals(arma_col1), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 0.004896##
## Box-Ljung test
##
## data: residuals(arma_col1)
## X-squared = 14.908, df = 4, p-value = 0.004896#ARMA (1,1)
arma_col11= Arima(d1col, c(1,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_col11) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.26782542 0.09846849 2.7199 0.00653 **
## ma1 -0.99999954 0.03003215 -33.2976 < 2e-16 ***
## drift -0.00021179 0.00014742 -1.4367 0.15081
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_col11) ## Series: d1col
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.2678 -1.00 -2e-04
## s.e. 0.0985 0.03 1e-04
##
## sigma^2 = 0.0005456: log likelihood = 230.94
## AIC=-453.87 AICc=-453.45 BIC=-443.49
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001545817 0.0228871 0.01710292 -24.20714 87.70007 0.8222211
## ACF1
## Training set -0.006224227Acf(residuals(arma_col11)) 
Box.test(residuals(arma_col11), lag = 5, type = c("Ljung-Box"), fitdf = 2) #p-value = 0.198##
## Box-Ljung test
##
## data: residuals(arma_col11)
## X-squared = 4.6653, df = 3, p-value = 0.198Acf(residuals(arma_col11)^2)
Box.test(residuals(arma_col11)^2, lag = 5, type = c("Ljung-Box"), fitdf =2) #p-value = 3.168e-05##
## Box-Ljung test
##
## data: residuals(arma_col11)^2
## X-squared = 23.505, df = 3, p-value = 3.168e-05Прогноз:
forecast_col11= forecast(arma_col11, h=30)
plot(forecast_col11, col='red')
#SUP-F
supf_col= Fstats(d1col ~ 1, from = 0.1, to = NULL)
plot(supf_col , alpha = 0.01) + lines(breakpoints(supf_col ))
## integer(0)sup= breakpoints(supf_col )
sup$breakpoints #10## [1] 10sctest(supf_col , type = "supF") #p-value = 0.003341##
## supF test
##
## data: supf_col
## sup.F = 14.769, p-value = 0.003341#ARMA модель на периоде (1,10)
arma_col_10 = d1col[seq(1,10)]
sup_11 = Fstats(arma_col_10 ~ 1, from = 0.1, to = NULL)## Warning in Fstats(arma_col_10 ~ 1, from = 0.1, to = NULL): 'from' changed (was
## too small)## Warning in Fstats(arma_col_10 ~ 1, from = 0.1, to = NULL): 'to' changed (was too
## large)plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #7## [1] 7sctest(sup_11 , type = "supF") #p-value = 0.5356##
## supF test
##
## data: sup_11
## sup.F = 2.7907, p-value = 0.5356#ARMA модель на периоде (10,60)
arma_col_60 = d1col[seq(10,60)]
sup_11 = Fstats(arma_col_60 ~ 1, from = 0.1, to = NULL)
plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #16## [1] 16sctest(sup_11 , type = "supF") #p-value = 0.1852##
## supF test
##
## data: sup_11
## sup.F = 6.1121, p-value = 0.1852#ARMA модель на периоде (60,100)
arma_col_100 = d1col[seq(60,100)]
sup_11 = Fstats(arma_col_100 ~ 1, from = 0.1, to = NULL)
plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #18## [1] 18sctest(sup_11 , type = "supF") #p-value = 0.853##
## supF test
##
## data: sup_11
## sup.F = 2.0013, p-value = 0.853Куба
eacf(d1cub)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x o x o o o o o o o o o o o
## 5 x x x x o o o o o o o o o o
## 6 x o o x x o o o o o o o o o
## 7 x o o x x o o o o o o o o o#ARMA (0,1)
arma_cub0= Arima(d1cub, c(0,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_cub0) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -9.9999e-01 3.4491e-02 -28.9924 <2e-16 ***
## drift 4.5555e-05 2.4982e-04 0.1823 0.8553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_cub0) ## Series: d1cub
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## -1.0000 0e+00
## s.e. 0.0345 2e-04
##
## sigma^2 = 0.004445: log likelihood = 126.33
## AIC=-246.65 AICc=-246.4 BIC=-238.87
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0001912944 0.0656615 0.04351983 -25.7911 163.0244 0.7888482
## ACF1
## Training set 0.344662Acf(residuals(arma_cub0)) 
Box.test(residuals(arma_cub0), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 0.004179##
## Box-Ljung test
##
## data: residuals(arma_cub0)
## X-squared = 15.267, df = 4, p-value = 0.004179#ARMA (1,0)
arma_cub1= Arima(d1cub, c(1,1,0), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_cub1) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.32781278 0.09438098 -3.4733 0.0005141 ***
## drift -0.00035803 0.00542148 -0.0660 0.9473471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_cub1) ## Series: d1cub
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## -0.3278 -0.0004
## s.e. 0.0944 0.0054
##
## sigma^2 = 0.005208: log likelihood = 120.73
## AIC=-235.45 AICc=-235.2 BIC=-227.67
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -8.834712e-05 0.07107672 0.05060327 -146.4925 312.3391 0.917244
## ACF1
## Training set -0.05072924Acf(residuals(arma_cub1)) 
Box.test(residuals(arma_cub1), lag = 5, type = c("Ljung-Box"), fitdf = 1) #p-value = 0.0141##
## Box-Ljung test
##
## data: residuals(arma_cub1)
## X-squared = 12.483, df = 4, p-value = 0.0141#ARMA (1,1)
arma_cub11= Arima(d1cub, c(1,1,1), include.constant=TRUE, method = c("CSS-ML"))
coeftest(arma_cub11) ##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 3.5164e-01 9.4994e-02 3.7017 0.0002142 ***
## ma1 -9.9995e-01 2.7995e-02 -35.7195 < 2.2e-16 ***
## drift 3.0568e-05 3.4144e-04 0.0895 0.9286632
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_cub11) ## Series: d1cub
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.3516 -0.9999 0e+00
## s.e. 0.0950 0.0280 3e-04
##
## sigma^2 = 0.003969: log likelihood = 132.8
## AIC=-257.6 AICc=-257.18 BIC=-247.22
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0004524551 0.06173053 0.04299539 -59.20589 198.4461 0.7793421
## ACF1
## Training set -0.007832507Acf(residuals(arma_cub11)) 
Box.test(residuals(arma_cub11), lag = 5, type = c("Ljung-Box"), fitdf = 2) #p-value = 0.3592##
## Box-Ljung test
##
## data: residuals(arma_cub11)
## X-squared = 3.2183, df = 3, p-value = 0.3592Acf(residuals(arma_cub11)^2)
Box.test(residuals(arma_cub11)^2, lag = 5, type = c("Ljung-Box"), fitdf =2) #p-value = 0.178##
## Box-Ljung test
##
## data: residuals(arma_cub11)^2
## X-squared = 4.9167, df = 3, p-value = 0.178Прогноз:
forecast_cub11= forecast(arma_cub11, h=30)
plot(forecast_cub11, col='red')
#SUP-F
supf_cub= Fstats(d1cub ~ 1, from = 0.1, to = NULL)
plot(supf_cub , alpha = 0.01) + lines(breakpoints(supf_cub ))
## integer(0)sup= breakpoints(supf_cub )
sup$breakpoints #15## [1] 15sctest(supf_cub , type = "supF") #p-value = 0.5377##
## supF test
##
## data: supf_cub
## sup.F = 3.4331, p-value = 0.5377#ARMA модель на периоде (1,15)
arma_cub_15 = d1cub[seq(1,15)]
sup_11 = Fstats(arma_cub_15 ~ 1, from = 0.1, to = NULL)## Warning in Fstats(arma_cub_15 ~ 1, from = 0.1, to = NULL): 'from' changed (was
## too small)## Warning in Fstats(arma_cub_15 ~ 1, from = 0.1, to = NULL): 'to' changed (was too
## large)plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #12## [1] 12sctest(sup_11 , type = "supF") #p-value = 0.02083##
## supF test
##
## data: sup_11
## sup.F = 10.654, p-value = 0.02083#ARMA модель на периоде (15,60)
arma_cub_60 = d1cub[seq(15,60)]
sup_11 = Fstats(arma_cub_60 ~ 1, from = 0.1, to = NULL)
plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #5## [1] 5sctest(sup_11 , type = "supF") #p-value = 0.9632##
## supF test
##
## data: sup_11
## sup.F = 1.4804, p-value = 0.9632#ARMA модель на периоде (60,100)
arma_cub_100 = d1cub[seq(60,100)]
sup_11 = Fstats(arma_cub_100 ~ 1, from = 0.1, to = NULL)
plot(sup_11 , alpha = 0.01) + lines(breakpoints(sup_11 ))
## integer(0)breakpoints(sup_11 )$breakpoints #17## [1] 17sctest(sup_11 , type = "supF") #p-value = 0.03771##
## supF test
##
## data: sup_11
## sup.F = 9.6725, p-value = 0.03771Аргентина поясняется от Боливии; Колумбия пояняется от Кубы, Бразилии; Куба поясняется от Боливии; Бразилия поясняется от Колумбии, Кубы, Аргентины;
ccf(d1arg, d1bol , lag.max = 30, type = c("correlation"))
grangertest(d1arg, d1bol, order = 9) #0.2191 является## Granger causality test
##
## Model 1: d1bol ~ Lags(d1bol, 1:9) + Lags(d1arg, 1:9)
## Model 2: d1bol ~ Lags(d1bol, 1:9)
## Res.Df Df F Pr(>F)
## 1 72
## 2 81 -9 1.3674 0.2191grangertest(d1bol, d1arg, order = 28) #0.4441 не является## Granger causality test
##
## Model 1: d1arg ~ Lags(d1arg, 1:28) + Lags(d1bol, 1:28)
## Model 2: d1arg ~ Lags(d1arg, 1:28)
## Res.Df Df F Pr(>F)
## 1 15
## 2 43 -28 0.8894 0.6191ccf(d1arg, d1bra , lag.max = 30, type = c("correlation"))
grangertest(d1arg, d1bra, order = 1) #0.7348 не является## Granger causality test
##
## Model 1: d1bra ~ Lags(d1bra, 1:1) + Lags(d1arg, 1:1)
## Model 2: d1bra ~ Lags(d1bra, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 0.1154 0.7348grangertest(d1bra, d1arg, order = 16) #0.0098 ** является## Granger causality test
##
## Model 1: d1arg ~ Lags(d1arg, 1:16) + Lags(d1bra, 1:16)
## Model 2: d1arg ~ Lags(d1arg, 1:16)
## Res.Df Df F Pr(>F)
## 1 51
## 2 67 -16 2.3793 0.009829 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ccf(d1arg, d1col , lag.max = 30, type = c("correlation"))
grangertest(d1arg, d1col, order = 1) #0.4026 является## Granger causality test
##
## Model 1: d1col ~ Lags(d1col, 1:1) + Lags(d1arg, 1:1)
## Model 2: d1col ~ Lags(d1col, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 0.7067 0.4026grangertest(d1col, d1arg, order = 20) #0.3384 не является## Granger causality test
##
## Model 1: d1arg ~ Lags(d1arg, 1:20) + Lags(d1col, 1:20)
## Model 2: d1arg ~ Lags(d1arg, 1:20)
## Res.Df Df F Pr(>F)
## 1 39
## 2 59 -20 1.1572 0.3384ccf(d1arg, d1cub , lag.max = 30, type = c("correlation"))
grangertest(d1arg, d1cub, order = 1) #0.7933 не является## Granger causality test
##
## Model 1: d1cub ~ Lags(d1cub, 1:1) + Lags(d1arg, 1:1)
## Model 2: d1cub ~ Lags(d1cub, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 0.069 0.7933grangertest(d1cub, d1arg, order = 9) #0.7038 является## Granger causality test
##
## Model 1: d1arg ~ Lags(d1arg, 1:9) + Lags(d1cub, 1:9)
## Model 2: d1arg ~ Lags(d1arg, 1:9)
## Res.Df Df F Pr(>F)
## 1 72
## 2 81 -9 0.7033 0.7038ccf(d1bra, d1col , lag.max = 30, type = c("correlation"))
grangertest(d1bra, d1col, order = 23) #0.08859 . является## Granger causality test
##
## Model 1: d1col ~ Lags(d1col, 1:23) + Lags(d1bra, 1:23)
## Model 2: d1col ~ Lags(d1col, 1:23)
## Res.Df Df F Pr(>F)
## 1 30
## 2 53 -23 1.6888 0.08859 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(d1col, d1bra, order = 1) #0.1215 является## Granger causality test
##
## Model 1: d1bra ~ Lags(d1bra, 1:1) + Lags(d1col, 1:1)
## Model 2: d1bra ~ Lags(d1bra, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 2.4402 0.1215ccf(d1bra, d1bol , lag.max = 18, type = c("correlation"))
grangertest(d1bra, d1bol, order = 1) #0.6703 не является## Granger causality test
##
## Model 1: d1bol ~ Lags(d1bol, 1:1) + Lags(d1bra, 1:1)
## Model 2: d1bol ~ Lags(d1bol, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 0.1823 0.6703grangertest(d1bol, d1bra, order = 15) #0.5182 не является## Granger causality test
##
## Model 1: d1bra ~ Lags(d1bra, 1:15) + Lags(d1bol, 1:15)
## Model 2: d1bra ~ Lags(d1bra, 1:15)
## Res.Df Df F Pr(>F)
## 1 54
## 2 69 -15 0.9495 0.5182ccf(d1bra, d1cub , lag.max = 18, type = c("correlation"))
grangertest(d1bra, d1cub, order = 1) #0.007944 ** является## Granger causality test
##
## Model 1: d1cub ~ Lags(d1cub, 1:1) + Lags(d1bra, 1:1)
## Model 2: d1cub ~ Lags(d1cub, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 12.006 0.0007944 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(d1cub, d1bra, order = 17) #0.7336 не является## Granger causality test
##
## Model 1: d1bra ~ Lags(d1bra, 1:17) + Lags(d1cub, 1:17)
## Model 2: d1bra ~ Lags(d1bra, 1:17)
## Res.Df Df F Pr(>F)
## 1 48
## 2 65 -17 0.7532 0.7336ccf(d1bol, d1cub , lag.max = 30, type = c("correlation"))
grangertest(d1bol, d1cub, order = 25) #0.4658 не является## Granger causality test
##
## Model 1: d1cub ~ Lags(d1cub, 1:25) + Lags(d1bol, 1:25)
## Model 2: d1cub ~ Lags(d1cub, 1:25)
## Res.Df Df F Pr(>F)
## 1 24
## 2 49 -25 1.0369 0.4658grangertest(d1cub, d1bol, order = 6) #0.04143 * является## Granger causality test
##
## Model 1: d1bol ~ Lags(d1bol, 1:6) + Lags(d1cub, 1:6)
## Model 2: d1bol ~ Lags(d1bol, 1:6)
## Res.Df Df F Pr(>F)
## 1 81
## 2 87 -6 2.3099 0.04143 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ccf(d1bol, d1col , lag.max = 30, type = c("correlation"))
grangertest(d1bol, d1col, order = 1) #0.7814 не является## Granger causality test
##
## Model 1: d1col ~ Lags(d1col, 1:1) + Lags(d1bol, 1:1)
## Model 2: d1col ~ Lags(d1col, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 0.0774 0.7814grangertest(d1col, d1bol, order = 16) #0.2204 не является## Granger causality test
##
## Model 1: d1bol ~ Lags(d1bol, 1:16) + Lags(d1col, 1:16)
## Model 2: d1bol ~ Lags(d1bol, 1:16)
## Res.Df Df F Pr(>F)
## 1 51
## 2 67 -16 1.3225 0.2204ccf(d1col, d1cub , lag.max = 50, type = c("correlation"))
grangertest(d1col, d1cub, order = 4) #0.162 является## Granger causality test
##
## Model 1: d1cub ~ Lags(d1cub, 1:4) + Lags(d1col, 1:4)
## Model 2: d1cub ~ Lags(d1cub, 1:4)
## Res.Df Df F Pr(>F)
## 1 87
## 2 91 -4 1.6794 0.162grangertest(d1cub, d1col, order = 1) #0.4936 не является## Granger causality test
##
## Model 1: d1col ~ Lags(d1col, 1:1) + Lags(d1cub, 1:1)
## Model 2: d1col ~ Lags(d1col, 1:1)
## Res.Df Df F Pr(>F)
## 1 96
## 2 97 -1 0.4722 0.4936Формула: BRA ~ COL + CUB + ARG
Оценим параметры ARMA модели для CPI и проверим ARCH-эффект:
eacf(d1bra)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o omodel_bra <- Arima(brazil$log_dyn, c(1,1,1), include.constant=TRUE, method = c("CSS"))
summary(model_bra)## Series: brazil$log_dyn
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.5858 -0.3297 0.0446
## s.e. 0.1799 0.1978 0.0065
##
## sigma^2 = 0.001618: log likelihood = 180.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.992876e-05 0.03942404 0.0296187 0.000463644 0.1523785 0.5472847
## ACF1
## Training set -0.03452445Acf(residuals(model_bra))
Box.test(residuals(model_bra), lag = 5, type = c("Ljung-Box"), fitdf = 3) #p-value = 0.2131##
## Box-Ljung test
##
## data: residuals(model_bra)
## X-squared = 3.0922, df = 2, p-value = 0.2131forecast_bra= forecast(model_bra, h=30)
plot(forecast_bra)
model_bra$aic #NA## [1] NA#Проверим ARCH эффект. ARMA(1,1)+ARCH 1
library(fGarch)## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
##
## If needed attach them yourself in your R script by e.g.,
## require("timeSeries")library(rugarch)##
## Attaching package: 'rugarch'## The following object is masked from 'package:stats':
##
## sigmagarch_bra= ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1, 0)),
mean.model= list(armaOrder = c(1,1), include.mean = TRUE, archm = FALSE), distribution.model = "norm")
garch_bra_fit = ugarchfit(garch_bra, d1bra)
garch_bra_fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,0)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu 0.047260 0.006939 6.8110 0.000000
## ar1 0.659723 0.206179 3.1998 0.001375
## ma1 -0.407530 0.201247 -2.0250 0.042864
## omega 0.001586 0.000359 4.4169 0.000010
## alpha1 0.000000 0.267369 0.0000 1.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu 0.047260 0.007125 6.63312 0.000000
## ar1 0.659723 0.279715 2.35855 0.018346
## ma1 -0.407530 0.673002 -0.60554 0.544819
## omega 0.001586 0.001842 0.86130 0.389075
## alpha1 0.000000 1.223394 0.00000 1.000000
##
## LogLikelihood : 179.0962
##
## Information Criteria
## ------------------------------------
##
## Akaike -3.4819
## Bayes -3.3517
## Shibata -3.4866
## Hannan-Quinn -3.4292
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.047 0.8284
## Lag[2*(p+q)+(p+q)-1][5] 1.022 1.0000
## Lag[4*(p+q)+(p+q)-1][9] 2.323 0.9662
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.001723 0.9669
## Lag[2*(p+q)+(p+q)-1][2] 0.008355 0.9912
## Lag[4*(p+q)+(p+q)-1][5] 0.104034 0.9980
## d.o.f=1
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[2] 0.01274 0.500 2.000 0.9101
## ARCH Lag[4] 0.09706 1.397 1.611 0.9846
## ARCH Lag[6] 0.16574 2.222 1.500 0.9972
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 1.1204
## Individual Statistics:
## mu 0.26206
## ar1 0.19892
## ma1 0.13554
## omega 0.07493
## alpha1 0.39489
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.28 1.47 1.88
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 1.0581 0.2927
## Negative Sign Bias 0.6254 0.5332
## Positive Sign Bias 0.7009 0.4851
## Joint Effect 1.1541 0.7640
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 16.8 0.6034
## 2 30 24.2 0.7190
## 3 40 36.0 0.6075
## 4 50 40.0 0.8168
##
##
## Elapsed time : 0.09349799coef(garch_bra_fit)## mu ar1 ma1 omega alpha1
## 4.726044e-02 6.597233e-01 -4.075304e-01 1.586424e-03 4.348173e-11Acf(residuals(garch_bra_fit))
Acf(residuals(garch_bra_fit)^2)
Acf(residuals(garch_bra_fit, standardize="TRUE"))
Acf(residuals(garch_bra_fit, standardize="TRUE")^2)
### Тест Bai Perron:
d11= ts(d1bra, start=1)
d111= c(0,d11[1:length(d11)-1])
brpoi= breakpoints(d11 ~ d111)
summary(brpoi)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = d11 ~ d111)
##
## Breakpoints at observation number:
##
## m = 1 62
## m = 2 24 62
## m = 3 24 39 62
## m = 4 24 39 62 81
## m = 5 24 39 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 62
## m = 2 24 62
## m = 3 24 39 62
## m = 4 24 39 62 81
## m = 5 24 39 55 70 85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 0.1688 0.1538 0.1414 0.1284 0.1281 0.1311
## BIC -340.8290 -336.3302 -330.9086 -326.7066 -313.1407 -297.0385plot(brpoi)
plot(d11) 
breakdates(brpoi)## [1] NA#con <- confint(brpoi)
#breakdates(con)
#con
#lines(con)d1_log<-diff(brazil$log_dyn, differences=1)
d1_log_l1 <- c(0,d1_log[1:length(d1_log)-1])
d1_log_l2 <- c(0,0,d1_log[2:length(d1_log)-2])SUP-F тест:
fstat_log= Fstats(d1_log ~ d1_log_l1, from = 0.1, to = NULL)
plot(fstat_log, alpha = 0.01) + lines(breakpoints(fstat_log))
## integer(0)kk= breakpoints(fstat_log)
kk$breakpoints #62## [1] 62sctest(fstat_log, type = "supF") #p-value = 0.1475 ##
## supF test
##
## data: fstat_log
## sup.F = 9.3739, p-value = 0.1475Изменение структуры OLS-CUSUM:
data_new= data.frame(d1_log, d1_log_l1, d1_log_l2)
efp_st= efp(d1_log ~ d1_log_l1, type = "OLS-CUSUM", data = data_new)
plot(efp_st, alpha = 0.1, functional = NULL)
sctest(efp_st) #p-value = 0.0611 разрыва нет##
## OLS-based CUSUM test
##
## data: efp_st
## S0 = 1.3207, p-value = 0.0611RE:
re_st <- efp(d1_log ~ d1_log_l1, type = "RE", data = data_new)
plot(re_st, alpha = 0.1, functional = NULL)
sctest(re_st) #p-value = 0.114 разрыва нет##
## RE test (recursive estimates test)
##
## data: re_st
## RE = 1.3282, p-value = 0.114ME:
me_st <- efp(d1_log ~ d1_log_l1, type = "ME",h=0.5, data = data_new)
plot(me_st, alpha = 0.1, functional = NULL)
sctest(me_st) #p-value = 0.1209 разрыва нет##
## ME test (moving estimates test)
##
## data: me_st
## ME = 1.4717, p-value = 0.1209Посмотрим на частичные структурные разрывы:
#partial SB xlag1
data_2new= data.frame(d1_log, d1_log_l1)
gefp_st= gefp(d1_log ~ d1_log_l1, data = data_2new, parm = "d1_log_l1")
plot(gefp_st, alpha = 0.1)
sctest(gefp_st) #p-value = 0.2866##
## M-fluctuation test
##
## data: gefp_st
## f(efp) = 0.98486, p-value = 0.2866Константа:
gefp_st_in <- gefp(d1_log ~ d1_log_l1, data = data_2new, parm = 1)
plot(gefp_st_in, alpha = 0.1)
sctest(gefp_st_in) #p-value = 0.06083##
## M-fluctuation test
##
## data: gefp_st_in
## f(efp) = 1.3215, p-value = 0.06083У нас есть частичный разрыв в константе. Разрыв в 62. Построим прогноз после разрыва.
ИЗМЕНЕНИЯ
Оцениваем параметры ARMA модели после разрыва:
eacf(d1bra)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o omodel_bra_razr <- Arima(log(brazil$log_dyn[62:101]), c(1,1,1), include.constant=TRUE, method = c("CSS"))
summary(model_bra_razr)## Series: log(brazil$log_dyn[62:101])
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## -0.6182 0.6432 0.0013
## s.e. 0.1910 0.2014 0.0003
##
## sigma^2 = 3.347e-06: log likelihood = 191.56
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set 1.616946e-05 0.001735686 0.001227119 0.0005215707 0.04013601
## MASE ACF1
## Training set 0.618165 0.1530751Acf(residuals(model_bra_razr))
Box.test(residuals(model_bra_razr), lag = 5, type = c("Ljung-Box"), fitdf = 0) #p-value = 0.3875##
## Box-Ljung test
##
## data: residuals(model_bra_razr)
## X-squared = 5.2379, df = 5, p-value = 0.3875forecast_razr= forecast(model_bra_razr, h=30)
plot(forecast_razr, col='red')
plot(forecast_bra)
Forecast без разрыва:
mean_arma= c(brazil$log_dyn[62:101], forecast_bra$mean)
up_arma= c(brazil$log_dyn[62:101], forecast_bra$upper[,2])
lower_arma= c(brazil$log_dyn[62:101], forecast_bra$lower[,2])
plot(mean_arma, lty=2, lwd=1, col='green',type='l', main= 'Прогноз реального ВВП Бразилии без разрыва', ylab='log(ВВП*популярность)')+ lines(up_arma, col='red',lty=2, lwd=1) + lines(lower_arma, col='purple',lty=2, lwd=1)
## integer(0)Forecast с разрывом:
mean_arma_razr= c(brazil$log_dyn[62:101], forecast_razr$mean)
up_arma_razr= c(brazil$log_dyn[62:101], forecast_razr$upper[,2])
lower_arma_razr= c(brazil$log_dyn[62:101], forecast_razr$lower[,2])
plot(mean_arma_razr)+ lines(mean_arma_razr, col='black')+
lines(up_arma_razr, col='red')+
lines(lower_arma, col='green')
## integer(0)Совместим 2 графика:
plot(mean_arma_razr, col='pink',type='l', main= 'Прогноз реального ВВП Бразилии с разрывом и без', ylab='log(ВВП*популярность)')+
lines(up_arma_razr, col='black',lty=2, lwd=1) +
lines(lower_arma_razr, col='blue',lty=2, lwd=1)+
lines(up_arma, col='red',lty=2, lwd=1) +
lines(lower_arma, col='purple',lty=2, lwd=1)+
lines(mean_arma, lty=2, lwd=1, col='green',type='l')
## integer(0)Прогнозируем популяцию. Добавляем новую колонку с лог(популяция) для Бразилии.
plot(brazil$year, brazil$pop,type='l',xlab='Год',ylab='кол-во, тыс',
main = 'Рост населения в Бразилии')
brazil= brazil %>% mutate(log_pop=log(brazil$pop))
d1bra_pop= diff(brazil$log_pop, differences = 1)
#Проверим стационарность
pacf(brazil$log_pop)
ur.df(brazil$log_pop,type='drift',selectlags = 'Fixed', lags=1) %>% summary()##
## ###############################################
## # 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.0178501 -0.0010933 -0.0001963 0.0009662 0.0277625
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.018236 0.007584 2.405 0.0181 *
## z.lag.1 -0.001298 0.000610 -2.129 0.0358 *
## z.diff.lag 0.825707 0.057595 14.337 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003527 on 96 degrees of freedom
## Multiple R-squared: 0.7658, Adjusted R-squared: 0.7609
## F-statistic: 156.9 on 2 and 96 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -2.1287 4.3363
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81pacf(d1bra_pop)
ur.df(d1bra_pop, type='drift',selectlags = 'Fixed', lags=0) %>% summary()##
## ###############################################
## # 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.0190757 -0.0007813 0.0000649 0.0004479 0.0282623
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002267 0.001134 1.999 0.0484 *
## z.lag.1 -0.114701 0.051241 -2.238 0.0275 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.00359 on 97 degrees of freedom
## Multiple R-squared: 0.04912, Adjusted R-squared: 0.03932
## F-statistic: 5.011 on 1 and 97 DF, p-value: 0.02748
##
##
## Value of test-statistic is: -2.2385 2.5804
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86Оценим параметры ARMA модели для популяции Бразилии:
eacf(brazil$log_pop)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x x x x
## 1 x x x x x x x x x x x x x x
## 2 x o o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x o o o o o o o o o o o o o
## 5 x o o o o o o o o o o o o o
## 6 x o o o o o o o o o o o o o
## 7 x o x o o o o o o o o o o oarma_bra_pop= Arima(brazil$log_pop, c(2,0,1), include.constant=TRUE, method = c("CSS"))
coeftest(arma_bra_pop)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.9736791 0.0065185 302.783 < 2.2e-16 ***
## ar2 -0.9741792 0.0064480 -151.082 < 2.2e-16 ***
## ma1 -0.6852399 0.0639508 -10.715 < 2.2e-16 ***
## intercept 12.2034496 0.2401687 50.812 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bra_pop)## Series: brazil$log_pop
## ARIMA(2,0,1) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 mean
## 1.9737 -0.9742 -0.6852 12.2034
## s.e. 0.0065 0.0064 0.0640 0.2402
##
## sigma^2 = 8.681e-06: log likelihood = 446.26
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set -9.758147e-05 0.002887414 0.0009588693 -0.000862064 0.008703773
## MASE ACF1
## Training set 0.04599728 -0.03519672Acf(residuals(arma_bra_pop))
Box.test(residuals(arma_bra_pop), lag = 5, type = c("Ljung-Box"), fitdf = 3) #p-value = 0.848##
## Box-Ljung test
##
## data: residuals(arma_bra_pop)
## X-squared = 0.3297, df = 2, p-value = 0.848forecast_pop= forecast(arma_bra_pop, h=30)
plot(forecast_pop)
plot(d1bra_pop, type='l')
Bai-Perron:
d1_pop_ts= ts(d1bra_pop, start=1)
d1_pop_ts_l1= c(0,d1_pop_ts[1:length(d1_pop_ts)-1])
brpoi_pop= breakpoints(d1_pop_ts ~ d1_pop_ts_l1)
summary(brpoi_pop)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = d1_pop_ts ~ d1_pop_ts_l1)
##
## Breakpoints at observation number:
##
## m = 1 33
## m = 2 31 46
## m = 3 16 31 46
## m = 4 16 31 46 64
## m = 5 16 31 46 64 83
##
## Corresponding to breakdates:
##
## m = 1 33
## m = 2 31 46
## m = 3 16 31 46
## m = 4 16 31 46 64
## m = 5 16 31 46 64 83
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 1.570e-03 9.091e-04 4.579e-04 4.101e-04 4.086e-04 4.072e-04
## BIC -8.086e+02 -8.494e+02 -9.042e+02 -9.014e+02 -8.879e+02 -8.745e+02plot(brpoi_pop)
Рассмотрим чистое ВВП Бразилии. Добавим колонку лог(ВВП).
brazil = brazil %>% mutate(log_gdp= log(brazil$gdppc))
#Разности
d1bra_gdp= diff(brazil$log_gdp, differences = 1)
#Проверка стационарности
pacf(brazil$log_gdp)
ur.df(brazil$log_gdp, type='drift',lags=0, selectlags = 'Fixed') %>% summary()##
## ###############################################
## # 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.151581 -0.020205 0.005161 0.028627 0.080155
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.064178 0.040624 1.580 0.117
## z.lag.1 -0.004582 0.004853 -0.944 0.347
##
## Residual standard error: 0.04093 on 98 degrees of freedom
## Multiple R-squared: 0.009017, Adjusted R-squared: -0.001095
## F-statistic: 0.8917 on 1 and 98 DF, p-value: 0.3473
##
##
## Value of test-statistic is: -0.9443 20.637
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81pacf(d1bra_gdp)
ur.df(d1bra_gdp,type='drift',lags=1, selectlags = 'Fixed') %>% summary()##
## ###############################################
## # 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.154557 -0.019954 0.004621 0.028081 0.094588
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.017067 0.005162 3.306 0.00134 **
## z.lag.1 -0.703197 0.126099 -5.577 2.3e-07 ***
## z.diff.lag -0.127683 0.099725 -1.280 0.20354
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03955 on 95 degrees of freedom
## Multiple R-squared: 0.4147, Adjusted R-squared: 0.4024
## F-statistic: 33.66 on 2 and 95 DF, p-value: 8.889e-12
##
##
## Value of test-statistic is: -5.5765 15.591
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.51 -2.89 -2.58
## phi1 6.70 4.71 3.86SUP-F:
d1_bra_gdp_l1 <- c(0,d1bra_gdp[1:length(d1bra_gdp)-1])
#Sup-F test
stat_bra_gdp= Fstats(d1bra_gdp ~ d1_bra_gdp_l1, from = 0.1, to = NULL)
plot(stat_bra_gdp, alpha = 0.01)+ lines(breakpoints(stat_bra_gdp))
## integer(0)breakpoints(stat_bra_gdp)$breakpoints #62## [1] 62sctest(stat_bra_gdp, type = "supF") #p-value = 0.2453##
## supF test
##
## data: stat_bra_gdp
## sup.F = 4.3394, p-value = 0.7402Bai-Perron:
d1_gdp_ts= ts(d1bra_gdp, start=1)
d1_gdp_ts_l1= c(0,d1_gdp_ts[1:length(d1_gdp_ts)-1])
stat_gdp= breakpoints(d1_gdp_ts ~ d1_gdp_ts_l1)
summary(stat_gdp)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = d1_gdp_ts ~ d1_gdp_ts_l1)
##
## Breakpoints at observation number:
##
## m = 1 62
## m = 2 22 39
## m = 3 22 39 62
## m = 4 22 39 62 81
## m = 5 16 32 47 62 81
##
## Corresponding to breakdates:
##
## m = 1 62
## m = 2 22 39
## m = 3 22 39 62
## m = 4 22 39 62 81
## m = 5 16 32 47 62 81
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 0.1591 0.1522 0.1424 0.1300 0.1283 0.1299
## BIC -346.7342 -337.3397 -330.1942 -325.4523 -313.0082 -297.9546plot(stat_gdp)
breakdates(stat_gdp) #NA## [1] NAОценим ARMA после 62:
d1bra_gdp_62= diff(brazil$log_gdp[62:101], differences = 1)
#Проверка стационарности:
pacf(d1bra_gdp_62)
ur.df(d1bra_gdp_62, type='drift', lags=1, selectlags = 'Fixed') %>% summary()##
## ###############################################
## # 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.147434 -0.014790 0.005969 0.024484 0.054352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.011197 0.007054 1.587 0.12169
## z.lag.1 -0.669597 0.209521 -3.196 0.00301 **
## z.diff.lag -0.190189 0.153559 -1.239 0.22400
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03841 on 34 degrees of freedom
## Multiple R-squared: 0.4751, Adjusted R-squared: 0.4442
## F-statistic: 15.39 on 2 and 34 DF, p-value: 1.743e-05
##
##
## Value of test-statistic is: -3.1958 5.1238
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.58 -2.93 -2.60
## phi1 7.06 4.86 3.94#eacf(d1bra_gdp_62)
arma_bra_gdp_62 <- Arima(brazil$log_gdp[62:101], c(0,1,1), include.constant=TRUE, method = c("CSS"))
coeftest(arma_bra_gdp_62)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.0606991 0.1330722 0.4561 0.64829
## drift 0.0153454 0.0068157 2.2515 0.02436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(arma_bra_gdp_62)## Series: brazil$log_gdp[62:101]
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.0607 0.0153
## s.e. 0.1331 0.0068
##
## sigma^2 = 0.0017: log likelihood = 70.04
##
## Training set error measures:
## ME RMSE MAE MPE MAPE
## Training set -8.193714e-05 0.03965683 0.02794663 -0.002240704 0.3034786
## MASE ACF1
## Training set 0.8266628 0.01313868Acf(residuals(arma_bra_gdp_62))
Box.test(residuals(arma_bra_gdp_62), lag = 6, type = c("Ljung-Box"), fitdf = 3) #p-value = 0.2679##
## Box-Ljung test
##
## data: residuals(arma_bra_gdp_62)
## X-squared = 3.9411, df = 3, p-value = 0.2679forecast_gdp_62= forecast(arma_bra_gdp_62, h=30)
plot(forecast_gdp_62)
Bai-Perron:
d1_gdp_ts_62= ts(d1bra_gdp_62, start=1)
d1_gdp_62_l1 <- c(0,d1_gdp_ts_62[1:length(d1_gdp_ts_62)-1])
stat_gdp_62 <- breakpoints(d1_gdp_ts_62 ~ d1_gdp_62_l1)
summary(stat_gdp_62)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = d1_gdp_ts_62 ~ d1_gdp_62_l1)
##
## Breakpoints at observation number:
##
## m = 1 34
## m = 2 13 34
## m = 3 5 13 34
## m = 4 5 11 24 34
## m = 5 5 13 18 24 34
## m = 6 5 13 18 24 29 34
##
## Corresponding to breakdates:
##
## m = 1 34
## m = 2 13 34
## m = 3 5 13 34
## m = 4 5 11 24 34
## m = 5 5 13 18 24 34
## m = 6 5 13 18 24 29 34
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 0.06289 0.05354 0.04657 0.04297 0.04139 0.04010
## BIC -129.09866 -124.38630 -118.83733 -110.97925 -101.45204 -91.69989
##
## m 6
## RSS 0.03916
## BIC -81.62923plot(stat_gdp_62)
forecast_mean= forecast_gdp_62$mean+forecast_pop$mean## Warning in .cbind.ts(list(e1, e2), c(deparse(substitute(e1))[1L],
## deparse(substitute(e2))[1L]), : non-intersecting seriesforecast_up= forecast_gdp_62$upper[,2]+forecast_pop$upper[,2]## Warning in .cbind.ts(list(e1, e2), c(deparse(substitute(e1))[1L],
## deparse(substitute(e2))[1L]), : non-intersecting seriesforecast_lower= forecast_gdp_62$lower[,2]+forecast_pop$lower[,2]## Warning in .cbind.ts(list(e1, e2), c(deparse(substitute(e1))[1L],
## deparse(substitute(e2))[1L]), : non-intersecting seriesbra_mean= c(brazil$log_dyn[62:101], forecast_mean)
bra_up= c(brazil$log_dyn[62:101], forecast_up)
bra_low= c(brazil$log_dyn[62:101], forecast_lower)
plot(bra_mean, col='red',type='l')+ lines(bra_up, col='green',type='l')+ lines(bra_low, col='pink',type='l')
## integer(0)Проврим коинтеграцию:
data_coin= data.frame(argentina$log_dyn, brazil$log_dyn, bolivia$log_dyn, colombia$log_dyn, cuba$log_dyn)
ca_coin= ca.jo(data_coin, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_coin)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.764457e-01 2.941769e-01 9.562769e-02 8.302815e-02 3.516378e-02
## [6] 8.709996e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 4 | 3.54 7.52 9.24 12.97
## r <= 3 | 8.58 13.75 15.67 20.20
## r <= 2 | 9.95 19.77 22.00 26.81
## r <= 1 | 34.49 25.56 28.14 33.24
## r = 0 | 46.76 31.66 34.40 39.79
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 bolivia.log_dyn.l1
## argentina.log_dyn.l1 1.00000000 1.00000000 1.0000000
## brazil.log_dyn.l1 0.72597122 -0.00262104 1.0415073
## bolivia.log_dyn.l1 -0.04522073 -0.28362351 1.0276410
## colombia.log_dyn.l1 -1.86245111 -0.29226463 -2.4499497
## cuba.log_dyn.l1 0.40422038 -0.29973515 -0.3495161
## constant -4.80511663 -4.03675541 -5.8736597
## colombia.log_dyn.l1 cuba.log_dyn.l1 constant
## argentina.log_dyn.l1 1.0000000 1.000000 1.00000
## brazil.log_dyn.l1 -0.9664726 6.580518 -42.78407
## bolivia.log_dyn.l1 -0.9509271 2.145256 19.02562
## colombia.log_dyn.l1 0.1009525 -9.912808 30.71002
## cuba.log_dyn.l1 1.5417872 2.539689 15.47158
## constant -12.7464732 -45.712238 -316.51100
##
## Weights W:
## (This is the loading matrix)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 bolivia.log_dyn.l1
## argentina.log_dyn.d -0.0158717064 -0.13847653 -0.031589495
## brazil.log_dyn.d 0.0351446203 0.15938961 -0.027491442
## bolivia.log_dyn.d 0.0008696852 0.06901314 -0.035543747
## colombia.log_dyn.d 0.0453260957 -0.04263678 -0.004685869
## cuba.log_dyn.d -0.0140217544 0.17670606 0.043390174
## colombia.log_dyn.l1 cuba.log_dyn.l1 constant
## argentina.log_dyn.d -0.023297693 -0.0020858722 2.337903e-14
## brazil.log_dyn.d -0.001810704 -0.0017678394 -1.029074e-14
## bolivia.log_dyn.d -0.018746207 0.0013548127 2.959728e-15
## colombia.log_dyn.d -0.007389471 -0.0002320196 -5.937995e-15
## cuba.log_dyn.d -0.042010399 -0.0002487478 -1.936443e-14#N= 5, 2 из 5 переменных коинтегрированыВычисляем не коинтегрированные переменные:
Уберем Аргентину:
dat_find= data.frame(brazil$log_dyn, bolivia$log_dyn, colombia$log_dyn, cuba$log_dyn)
ca_find= ca.jo(dat_find, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_find)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.466261e-01 1.278774e-01 8.295275e-02 1.797030e-02 -3.947453e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 1.80 7.52 9.24 12.97
## r <= 2 | 8.57 13.75 15.67 20.20
## r <= 1 | 13.55 19.77 22.00 26.81
## r = 0 | 42.13 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## brazil.log_dyn.l1 bolivia.log_dyn.l1 colombia.log_dyn.l1
## brazil.log_dyn.l1 1.0000000 1.000000 1.0000000
## bolivia.log_dyn.l1 0.3706559 1.609539 0.9105622
## colombia.log_dyn.l1 -2.1140654 -3.038563 -1.1949314
## cuba.log_dyn.l1 0.8269432 1.061020 -0.9474777
## constant -0.1407040 -8.838856 3.3670935
## cuba.log_dyn.l1 constant
## brazil.log_dyn.l1 1.0000000 1.0000000
## bolivia.log_dyn.l1 0.1503642 -0.4481402
## colombia.log_dyn.l1 -1.6276119 -0.7367396
## cuba.log_dyn.l1 1.3140544 -0.3448023
## constant -14.6237274 7.0677120
##
## Weights W:
## (This is the loading matrix)
##
## brazil.log_dyn.l1 bolivia.log_dyn.l1 colombia.log_dyn.l1
## brazil.log_dyn.d 0.018307444 -0.026581402 -0.015424831
## bolivia.log_dyn.d -0.003823635 -0.025946104 0.010777078
## colombia.log_dyn.d 0.032563996 -0.004064995 0.005069904
## cuba.log_dyn.d -0.013873792 -0.008142669 0.087686623
## cuba.log_dyn.l1 constant
## brazil.log_dyn.d -0.0043628385 6.392687e-13
## bolivia.log_dyn.d 0.0041747814 5.239751e-13
## colombia.log_dyn.d 0.0002977997 -4.371996e-13
## cuba.log_dyn.d -0.0035871512 -4.834370e-13#1 из 4 коинтегрированыУберем Бразилию:
dat_find= data.frame(argentina$log_dyn, bolivia$log_dyn, colombia$log_dyn, cuba$log_dyn)
ca_find= ca.jo(dat_find, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_find)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.711908e-01 1.756004e-01 8.729999e-02 4.660482e-02 -3.199501e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 4.72 7.52 9.24 12.97
## r <= 2 | 9.04 13.75 15.67 20.20
## r <= 1 | 19.12 19.77 22.00 26.81
## r = 0 | 45.93 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## argentina.log_dyn.l1 bolivia.log_dyn.l1
## argentina.log_dyn.l1 1.0000000 1.00000000
## bolivia.log_dyn.l1 -0.1446572 -0.08227612
## colombia.log_dyn.l1 -1.6081095 -0.38512469
## cuba.log_dyn.l1 1.4191043 -0.35496236
## constant -10.2774005 -4.75625084
## colombia.log_dyn.l1 cuba.log_dyn.l1 constant
## argentina.log_dyn.l1 1.0000000 1.0000000 1.00000000
## bolivia.log_dyn.l1 0.2644479 -1.4770237 -2.40791294
## colombia.log_dyn.l1 -1.2046157 -0.3838997 0.79383480
## cuba.log_dyn.l1 0.4778794 1.6311003 -0.06587414
## constant -9.7214816 -15.6506043 6.96701391
##
## Weights W:
## (This is the loading matrix)
##
## argentina.log_dyn.l1 bolivia.log_dyn.l1 colombia.log_dyn.l1
## argentina.log_dyn.d -0.0044746136 -0.096165766 -0.05311398
## bolivia.log_dyn.d -0.0003030808 0.051025523 -0.05390050
## colombia.log_dyn.d 0.0265886057 -0.001576396 -0.01170599
## cuba.log_dyn.d -0.0064077600 0.255469274 -0.01176914
## cuba.log_dyn.l1 constant
## argentina.log_dyn.d -0.011350476 1.064888e-13
## bolivia.log_dyn.d 0.006213695 1.171717e-14
## colombia.log_dyn.d -0.002345878 -1.719558e-14
## cuba.log_dyn.d -0.011802375 -9.828556e-14#1 из 4 коинтегрированыУберем Боливию:
dat_find= data.frame(argentina$log_dyn, brazil$log_dyn, colombia$log_dyn, cuba$log_dyn)
ca_find= ca.jo(dat_find, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_find)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.791665e-01 2.407299e-01 7.015145e-02 1.937275e-02 -7.693614e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 1.94 7.52 9.24 12.97
## r <= 2 | 7.20 13.75 15.67 20.20
## r <= 1 | 27.26 19.77 22.00 26.81
## r = 0 | 47.19 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 colombia.log_dyn.l1
## argentina.log_dyn.l1 1.0000000 1.0000000 1.00000000
## brazil.log_dyn.l1 0.7449010 0.1328680 -0.04443937
## colombia.log_dyn.l1 -1.8701857 -0.5956987 -1.16865583
## cuba.log_dyn.l1 0.3462424 -0.3793071 0.88622154
## constant -4.7877018 -4.4718113 -12.09369528
## cuba.log_dyn.l1 constant
## argentina.log_dyn.l1 1.000000 1.000000
## brazil.log_dyn.l1 -14.860127 2.576792
## colombia.log_dyn.l1 13.900075 -3.601983
## cuba.log_dyn.l1 5.094968 -1.283543
## constant -68.163149 17.854411
##
## Weights W:
## (This is the loading matrix)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 colombia.log_dyn.l1
## argentina.log_dyn.d -0.014868108 -0.12911885 -0.05144538
## brazil.log_dyn.d 0.034900061 0.09899516 -0.01979231
## colombia.log_dyn.d 0.045776726 -0.03599840 -0.01078980
## cuba.log_dyn.d -0.009807937 0.23808751 -0.03849926
## cuba.log_dyn.l1 constant
## argentina.log_dyn.d 9.714594e-04 1.580994e-15
## brazil.log_dyn.d 2.603137e-03 -2.431387e-15
## colombia.log_dyn.d -9.799027e-05 -5.432592e-15
## cuba.log_dyn.d -2.721887e-03 -7.686405e-15#2 из 4 коинтегрированыУберем Колумбию:
dat_find= data.frame(argentina$log_dyn, brazil$log_dyn, bolivia$log_dyn, cuba$log_dyn)
ca_find= ca.jo(dat_find, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_find)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.288181e-01 1.488584e-01 8.802556e-02 4.272843e-02 -1.406182e-15
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 4.32 7.52 9.24 12.97
## r <= 2 | 9.12 13.75 15.67 20.20
## r <= 1 | 15.96 19.77 22.00 26.81
## r = 0 | 39.47 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 bolivia.log_dyn.l1
## argentina.log_dyn.l1 1.0000000 1.0000000 1.0000000
## brazil.log_dyn.l1 -0.2043660 -0.2558021 -0.8784797
## bolivia.log_dyn.l1 -0.5649930 -0.1889332 -0.7553958
## cuba.log_dyn.l1 -0.1060961 -0.3422399 1.3469399
## constant -3.8665406 -5.3230966 -12.4687814
## cuba.log_dyn.l1 constant
## argentina.log_dyn.l1 1.0000000 1.00000
## brazil.log_dyn.l1 -1.1106715 -20.52121
## bolivia.log_dyn.l1 0.9893285 26.88579
## cuba.log_dyn.l1 -0.6704074 13.35746
## constant -2.2726931 -284.10616
##
## Weights W:
## (This is the loading matrix)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 bolivia.log_dyn.l1
## argentina.log_dyn.d 0.010869368 -0.125525584 -0.036361821
## brazil.log_dyn.d 0.162362810 -0.031894611 -0.004461834
## bolivia.log_dyn.d 0.068045621 -0.007028568 -0.024865976
## cuba.log_dyn.d -0.004378703 0.132477154 -0.047053384
## cuba.log_dyn.l1 constant
## argentina.log_dyn.d 0.004011598 -4.956631e-15
## brazil.log_dyn.d 0.004451809 3.174753e-15
## bolivia.log_dyn.d -0.008749622 7.978354e-16
## cuba.log_dyn.d 0.005887784 5.429956e-15#1 из 4 коинтегрированыУберем Кубу:
dat_find= data.frame(argentina$log_dyn, brazil$log_dyn, colombia$log_dyn, bolivia$log_dyn)
ca_find= ca.jo(dat_find, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_find)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 3.608256e-01 2.331903e-01 7.968863e-02 3.817126e-02 1.288403e-15
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.85 7.52 9.24 12.97
## r <= 2 | 8.22 13.75 15.67 20.20
## r <= 1 | 26.29 19.77 22.00 26.81
## r = 0 | 44.31 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 colombia.log_dyn.l1
## argentina.log_dyn.l1 1.0000000 1.00000000 1.0000000
## brazil.log_dyn.l1 0.4600608 -0.34064359 0.8646964
## colombia.log_dyn.l1 -1.1512225 0.03253577 -2.1503946
## bolivia.log_dyn.l1 -0.1921143 -0.42673295 0.5447174
## constant -3.6093532 -6.12593623 -5.8875575
## bolivia.log_dyn.l1 constant
## argentina.log_dyn.l1 1.00000 1.000000
## brazil.log_dyn.l1 20.75088 -4.670597
## colombia.log_dyn.l1 -29.23472 2.082178
## bolivia.log_dyn.l1 12.32450 4.685350
## constant -94.89602 -42.589855
##
## Weights W:
## (This is the loading matrix)
##
## argentina.log_dyn.l1 brazil.log_dyn.l1 colombia.log_dyn.l1
## argentina.log_dyn.d -0.03887569 -0.12218045 -0.047312174
## brazil.log_dyn.d 0.07295043 0.09557829 -0.028346898
## colombia.log_dyn.d 0.06982500 -0.06744200 -0.008159958
## bolivia.log_dyn.d 0.01236250 0.03575203 -0.044049288
## bolivia.log_dyn.l1 constant
## argentina.log_dyn.d -7.075316e-04 1.621915e-13
## brazil.log_dyn.d -6.971502e-04 -8.424571e-14
## colombia.log_dyn.d -7.545225e-05 -1.974531e-14
## bolivia.log_dyn.d 6.146090e-04 9.784474e-15#1 из 4 коинтегрированыУберем Боливию, Кубу, Колумбию. Для них построим VAR-модели.
dat_find= data.frame(brazil$log_dyn, argentina$log_dyn)
ca_find= ca.jo(dat_find, ecdet = "const", type="eigen", K=2, spec="transitory", season=NULL)
summary(ca_find)##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 2.456041e-01 8.666141e-02 9.893211e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 1 | 8.97 7.52 9.24 12.97
## r = 0 | 27.90 13.75 15.67 20.20
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## brazil.log_dyn.l1 argentina.log_dyn.l1 constant
## brazil.log_dyn.l1 1.000000 1.000000 1.0000000
## argentina.log_dyn.l1 -1.655988 -1.859185 -0.9775630
## constant 11.157339 16.118021 -0.7543483
##
## Weights W:
## (This is the loading matrix)
##
## brazil.log_dyn.l1 argentina.log_dyn.l1 constant
## brazil.log_dyn.d -0.03379174 -0.0004687576 -2.032305e-15
## argentina.log_dyn.d -0.01240336 0.0690107230 1.468738e-15#N = 4 (2 из 4 коинтегрированы)library(tsDyn)
vecm_find= VECM(dat_find, lag=2, include = c("none"), estim="ML", LRinclude = c("const"))
summary(vecm_find)## #############
## ###Model VECM
## #############
## Full sample size: 101 End sample size: 98
## Number of variables: 2 Number of estimated slope parameters 10
## AIC -1217.534 BIC -1189.099 SSR 0.375977
## Cointegrating vector (estimated by ML):
## brazil.log_dyn argentina.log_dyn const
## r1 1 -1.463032 6.366698
##
##
## ECT brazil.log_dyn -1
## Equation brazil.log_dyn -0.0136(0.0034)*** 0.2006(0.1044).
## Equation argentina.log_dyn -0.0111(0.0042)** 0.2616(0.1273)*
## argentina.log_dyn -1 brazil.log_dyn -2
## Equation brazil.log_dyn -0.0044(0.0875) 0.1596(0.1041)
## Equation argentina.log_dyn 0.0974(0.1067) -0.1350(0.1270)
## argentina.log_dyn -2
## Equation brazil.log_dyn 0.0213(0.0869)
## Equation argentina.log_dyn -0.1135(0.1060)vecm_find_pred = predict(vecm_find, n.ahead=30)
plot(vecm_find_pred)
prog1= c(brazil$log_dyn, vecm_find_pred[,1])
plot(prog1, type = 'l')+ title('Прогноз ВВП Бразилии на 30 лет')
## integer(0)vecm_find_pred= c(argentina$log_dyn, vecm_find_pred[,2])
plot(vecm_find_pred, type = 'l')+title('Прогноз ВВП Аргентины на 30 лет')
## integer(0)Строим VAR модель:
Бразилия и Колумбия:
d1_bra_col = data.frame(d1bra, d1col)
kk2<-VARselect(d1_bra_col, lag.max = 24, type="const")
var_mo= VAR(d1_bra_col, p = 2, type = "const")
summary(var_mo)##
## VAR Estimation Results:
## =========================
## Endogenous variables: d1bra, d1col
## Deterministic variables: const
## Sample size: 98
## Log Likelihood: 416.641
## Roots of the characteristic polynomial:
## 0.5105 0.4368 0.2408 0.05894
## Call:
## VAR(y = d1_bra_col, p = 2, type = "const")
##
##
## Estimation results for equation d1bra:
## ======================================
## d1bra = d1bra.l1 + d1col.l1 + d1bra.l2 + d1col.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.188231 0.107240 1.755 0.0825 .
## d1col.l1 0.229905 0.190563 1.206 0.2307
## d1bra.l2 0.191820 0.108502 1.768 0.0804 .
## d1col.l2 -0.193498 0.186468 -1.038 0.3021
## const 0.026321 0.009894 2.660 0.0092 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.04027 on 93 degrees of freedom
## Multiple R-Squared: 0.1244, Adjusted R-squared: 0.08671
## F-statistic: 3.302 on 4 and 93 DF, p-value: 0.0141
##
##
## Estimation results for equation d1col:
## ======================================
## d1col = d1bra.l1 + d1col.l1 + d1bra.l2 + d1col.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.144407 0.061387 2.352 0.0208 *
## d1col.l1 0.185194 0.109084 1.698 0.0929 .
## d1bra.l2 -0.012948 0.062109 -0.208 0.8353
## d1col.l2 -0.003438 0.106740 -0.032 0.9744
## const 0.028214 0.005664 4.982 2.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.02305 on 93 degrees of freedom
## Multiple R-Squared: 0.1291, Adjusted R-squared: 0.09164
## F-statistic: 3.447 on 4 and 93 DF, p-value: 0.01132
##
##
##
## Covariance matrix of residuals:
## d1bra d1col
## d1bra 0.0016215 0.0002987
## d1col 0.0002987 0.0005313
##
## Correlation matrix of residuals:
## d1bra d1col
## d1bra 1.0000 0.3218
## d1col 0.3218 1.0000coef(var_mo)## $d1bra
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.18823076 0.107240436 1.755222 0.082514346
## d1col.l1 0.22990481 0.190563121 1.206450 0.230704203
## d1bra.l2 0.19182009 0.108501854 1.767897 0.080357752
## d1col.l2 -0.19349845 0.186468219 -1.037702 0.302099350
## const 0.02632101 0.009894301 2.660220 0.009197364
##
## $d1col
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.144407228 0.061387406 2.3523918 2.076176e-02
## d1col.l1 0.185194066 0.109083626 1.6977256 9.290351e-02
## d1bra.l2 -0.012948027 0.062109476 -0.2084710 8.353171e-01
## d1col.l2 -0.003437559 0.106739590 -0.0322051 9.743775e-01
## const 0.028214378 0.005663773 4.9815516 2.900827e-06Acf(var_mo$varresult$d1bra$residuals)
Acf(var_mo$varresult$d1col$residuals)
ccf(var_mo$varresult$d1bra$residuals, var_mo$varresult$d1col$residuals)
#### Проведем тесты Hosking и LiMcLeod:
Hosking(var_mo, lags=2*var_mo$p)## lags statistic df p-value
## 4 10.81612 8 0.2123375LiMcLeod(var_mo, lags=2*var_mo$p)## lags statistic df p-value
## 4 10.87381 8 0.2089523Модель ARIMA:
eacf(d1bra)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o omodel_bra2 <- Arima(brazil$log_dyn, c(1,1,1), include.constant=TRUE, method = c("CSS"))
summary(model_bra2)## Series: brazil$log_dyn
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.5858 -0.3297 0.0446
## s.e. 0.1799 0.1978 0.0065
##
## sigma^2 = 0.001618: log likelihood = 180.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.992876e-05 0.03942404 0.0296187 0.000463644 0.1523785 0.5472847
## ACF1
## Training set -0.03452445Acf(residuals(model_bra2))
Box.test(residuals(model_bra2), lag = 5, type = c("Ljung-Box"), fitdf = length(coef(model_bra2))-1)##
## Box-Ljung test
##
## data: residuals(model_bra2)
## X-squared = 3.0922, df = 3, p-value = 0.3776Модель без ограничений:
ess= sum(var_mo$varresult$d1bra$residuals^2)
n1=var_mo$varresult$d1bra$rank-1
n= length(coef(model_bra2))-1
x= var_mo$p-n+1
x1= length(d1bra)-2*n1-1
rss= sum(model_bra2$residuals[x:length(residuals(model_bra2))]^2)
Fstat= ((ess-rss)/(n-n1))/(rss/x1)
Fstat #1.792431## [1] 1.792431qf(0.99, df1=n1-n, df2=length(d1bra)-2*n1-1) #4.846286## [1] 4.846286pf(Fstat, n1-n, length(d1bra)-2*n1-1, lower.tail=F) #0.1723848 ## [1] 0.1723848Бразилия и Куба
d1_bra_cub = data.frame(d1bra, d1cub)
kk2<-VARselect(d1_bra_cub, lag.max = 24, type="const")
var_mo= VAR(d1_bra_cub, p = 2, type = "const")
summary(var_mo)##
## VAR Estimation Results:
## =========================
## Endogenous variables: d1bra, d1cub
## Deterministic variables: const
## Sample size: 98
## Log Likelihood: 316.833
## Roots of the characteristic polynomial:
## 0.5131 0.3809 0.3445 0.03815
## Call:
## VAR(y = d1_bra_cub, p = 2, type = "const")
##
##
## Estimation results for equation d1bra:
## ======================================
## d1bra = d1bra.l1 + d1cub.l1 + d1bra.l2 + d1cub.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.216299 0.101363 2.134 0.035481 *
## d1cub.l1 -0.068463 0.069784 -0.981 0.329102
## d1bra.l2 0.216931 0.106107 2.044 0.043734 *
## d1cub.l2 0.013961 0.065718 0.212 0.832227
## const 0.026655 0.007128 3.740 0.000318 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.04054 on 93 degrees of freedom
## Multiple R-Squared: 0.1126, Adjusted R-squared: 0.07447
## F-statistic: 2.951 on 4 and 93 DF, p-value: 0.02409
##
##
## Estimation results for equation d1cub:
## ======================================
## d1cub = d1bra.l1 + d1cub.l1 + d1bra.l2 + d1cub.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.498606 0.150118 3.321 0.00128 **
## d1cub.l1 0.295113 0.103350 2.855 0.00530 **
## d1bra.l2 0.041431 0.157144 0.264 0.79263
## d1cub.l2 0.014508 0.097329 0.149 0.88183
## const -0.007088 0.010556 -0.671 0.50361
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.06003 on 93 degrees of freedom
## Multiple R-Squared: 0.2206, Adjusted R-squared: 0.1871
## F-statistic: 6.581 on 4 and 93 DF, p-value: 0.0001043
##
##
##
## Covariance matrix of residuals:
## d1bra d1cub
## d1bra 1.643e-03 3.454e-05
## d1cub 3.454e-05 3.604e-03
##
## Correlation matrix of residuals:
## d1bra d1cub
## d1bra 1.00000 0.01419
## d1cub 0.01419 1.00000coef(var_mo)## $d1bra
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.21629893 0.101362713 2.1339102 0.0354812725
## d1cub.l1 -0.06846310 0.069783876 -0.9810733 0.3291019630
## d1bra.l2 0.21693146 0.106106579 2.0444676 0.0437341101
## d1cub.l2 0.01396132 0.065718271 0.2124420 0.8322271626
## const 0.02665490 0.007127712 3.7396158 0.0003183845
##
## $d1cub
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.498605506 0.15011848 3.3214133 0.00128104
## d1cub.l1 0.295113036 0.10335013 2.8554685 0.00529959
## d1bra.l2 0.041430929 0.15714416 0.2636492 0.79263356
## d1cub.l2 0.014507884 0.09732895 0.1490603 0.88182881
## const -0.007087722 0.01055616 -0.6714298 0.50361115Acf(var_mo$varresult$d1bra$residuals)
Acf(var_mo$varresult$d1cub$residuals)
ccf(var_mo$varresult$d1bra$residuals, var_mo$varresult$d1cub$residuals)
#### Проведем тесты Hosking и LiMcLeod:
Hosking(var_mo, lags=2*var_mo$p)## lags statistic df p-value
## 4 3.987882 8 0.8582151LiMcLeod(var_mo, lags=2*var_mo$p)## lags statistic df p-value
## 4 4.258444 8 0.8330854Модель ARIMA:
eacf(d1bra)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o omodel_bra2 <- Arima(brazil$log_dyn, c(1,1,1), include.constant=TRUE, method = c("CSS"))
summary(model_bra2)## Series: brazil$log_dyn
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.5858 -0.3297 0.0446
## s.e. 0.1799 0.1978 0.0065
##
## sigma^2 = 0.001618: log likelihood = 180.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.992876e-05 0.03942404 0.0296187 0.000463644 0.1523785 0.5472847
## ACF1
## Training set -0.03452445Acf(residuals(model_bra2))
Box.test(residuals(model_bra2), lag = 12, type = c("Ljung-Box"), fitdf = length(coef(model_bra2))-1)##
## Box-Ljung test
##
## data: residuals(model_bra2)
## X-squared = 6.8878, df = 10, p-value = 0.736Модель без ограничений:
ess= sum(var_mo$varresult$d1bra$residuals^2)
n1=var_mo$varresult$d1bra$rank-1
n= length(coef(model_bra2))-1
x= var_mo$p-n+1
x1= length(d1bra)-2*n1-1
rss= sum(model_bra2$residuals[x:length(residuals(model_bra2))]^2)
Fstat= ((ess-rss)/(n-n1))/(rss/x1)
Fstat #1.20659## [1] 1.20659qf(0.99, df1=n1-n, df2=length(d1bra)-2*n1-1) #4.846286## [1] 4.846286pf(Fstat, n1-n, length(d1bra)-2*n1-1, lower.tail=F) #0.3039571## [1] 0.3039571Бразилия и Боливия:
d1_bra_bol = data.frame(d1bra, d1bol)
kk2<-VARselect(d1_bra_bol, lag.max = 24, type="const")
var_mo= VAR(d1_bra_bol, p = 2, type = "const")
summary(var_mo)##
## VAR Estimation Results:
## =========================
## Endogenous variables: d1bra, d1bol
## Deterministic variables: const
## Sample size: 98
## Log Likelihood: 363.339
## Roots of the characteristic polynomial:
## 0.4762 0.4762 0.3264 0.09111
## Call:
## VAR(y = d1_bra_bol, p = 2, type = "const")
##
##
## Estimation results for equation d1bra:
## ======================================
## d1bra = d1bra.l1 + d1bol.l1 + d1bra.l2 + d1bol.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.204063 0.101833 2.004 0.0480 *
## d1bol.l1 -0.003984 0.108760 -0.037 0.9709
## d1bra.l2 0.220433 0.100189 2.200 0.0303 *
## d1bol.l2 -0.205342 0.108630 -1.890 0.0618 .
## const 0.031713 0.007411 4.279 4.56e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.0399 on 93 degrees of freedom
## Multiple R-Squared: 0.1403, Adjusted R-squared: 0.1034
## F-statistic: 3.796 on 4 and 93 DF, p-value: 0.006647
##
##
## Estimation results for equation d1bol:
## ======================================
## d1bol = d1bra.l1 + d1bol.l1 + d1bra.l2 + d1bol.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.030482 0.098585 0.309 0.75787
## d1bol.l1 0.302420 0.105291 2.872 0.00505 **
## d1bra.l2 0.064317 0.096994 0.663 0.50890
## d1bol.l2 -0.029322 0.105166 -0.279 0.78100
## const 0.017198 0.007175 2.397 0.01853 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.03863 on 93 degrees of freedom
## Multiple R-Squared: 0.1012, Adjusted R-squared: 0.06258
## F-statistic: 2.619 on 4 and 93 DF, p-value: 0.0399
##
##
##
## Covariance matrix of residuals:
## d1bra d1bol
## d1bra 0.0015919 0.0002881
## d1bol 0.0002881 0.0014920
##
## Correlation matrix of residuals:
## d1bra d1bol
## d1bra 1.000 0.187
## d1bol 0.187 1.000coef(var_mo)## $d1bra
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.20406292 0.101832624 2.00390513 4.798951e-02
## d1bol.l1 -0.00398354 0.108759598 -0.03662702 9.708609e-01
## d1bra.l2 0.22043270 0.100189482 2.20015813 3.027289e-02
## d1bol.l2 -0.20534186 0.108630399 -1.89027988 6.183448e-02
## const 0.03171254 0.007411449 4.27885875 4.555834e-05
##
## $d1bol
## Estimate Std. Error t value Pr(>|t|)
## d1bra.l1 0.03048162 0.098584883 0.3091917 0.757867386
## d1bol.l1 0.30241990 0.105290935 2.8722311 0.005048383
## d1bra.l2 0.06431741 0.096994146 0.6631061 0.508902876
## d1bol.l2 -0.02932247 0.105165857 -0.2788212 0.781001093
## const 0.01719814 0.007175076 2.3969273 0.018533105Acf(var_mo$varresult$d1bra$residuals)
Acf(var_mo$varresult$d1bol$residuals)
ccf(var_mo$varresult$d1bra$residuals, var_mo$varresult$d1bol$residuals)
#### Проведем тесты Hosking и LiMcLeod:
Hosking(var_mo, lags=2*var_mo$p)## lags statistic df p-value
## 4 24.62764 8 0.001797146LiMcLeod(var_mo, lags=2*var_mo$p)## lags statistic df p-value
## 4 24.30702 8 0.002035222Модель ARIMA:
eacf(d1bra)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x o o o o o o o o o o o o
## 5 x x o o o o o o o o o o o o
## 6 x x o o o o o o o o o o o o
## 7 x x o o o o o o o o o o o omodel_bra2 <- Arima(brazil$log_dyn, c(1,1,1), include.constant=TRUE, method = c("CSS"))
summary(model_bra2)## Series: brazil$log_dyn
## ARIMA(1,1,1) with drift
##
## Coefficients:
## ar1 ma1 drift
## 0.5858 -0.3297 0.0446
## s.e. 0.1799 0.1978 0.0065
##
## sigma^2 = 0.001618: log likelihood = 180.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -3.992876e-05 0.03942404 0.0296187 0.000463644 0.1523785 0.5472847
## ACF1
## Training set -0.03452445Acf(residuals(model_bra2))
Box.test(residuals(model_bra2), lag = 5, type = c("Ljung-Box"), fitdf = length(coef(model_bra2))-1)##
## Box-Ljung test
##
## data: residuals(model_bra2)
## X-squared = 3.0922, df = 3, p-value = 0.3776Модель без ограничений:
ess= sum(var_mo$varresult$d1bra$residuals^2)
n1=var_mo$varresult$d1bra$rank-1
n= length(coef(model_bra2))-1
x= var_mo$p-n+1
x1= length(d1bra)-2*n1-1
rss= sum(model_bra2$residuals[x:length(residuals(model_bra2))]^2)
Fstat= ((ess-rss)/(n-n1))/(rss/x1)
Fstat #2.589876## [1] 2.589876qf(0.99, df1=n1-n, df2=length(d1bra)-2*n1-1) #4.846286## [1] 4.846286pf(Fstat, n1-n, length(d1bra)-2*n1-1, lower.tail=F) #0.08055217## [1] 0.08055217warnings()