Appendix for Natural Disaster Risk
RStudio Econometric Code
Nicolas Boccard
24 of April 2020
You can run this code directly at RStudio Cloud after a free registration and importing it into your own space. On the right panel named FILES, click on NDRiskCoding.Rmd to open it. Now you can tweak the code and run it by yourself but be patient the first time since a number of components will need to be installed. The Excel source is included.
1 Preliminaries
As a general rule, you can execute chunks of R code by clicking the Run button within the chunk. The first time you open this file, you’ll need to install the following packages (remove the comment sign # before executing). If you plan on using this RMD file on your PC, remove the # in front of the working directory address and update it to the correct local path.
# install.packages("tseries")
# install.packages("forecast")
# install.packages("FitAR")
# install.packages("strucchange")
# install.packages("readxl")
#setwd("~/Dropbox/Doks/Risque/NDtrend/data/")First, we load libraries and call the dataset, declaring its time series nature. To test later for trend breaks, we shall use a function developped for the Bai-Perron test implementation at https://rdrr.io/cran/strucchange/man/RealInt.html.
library(tseries)
library(forecast)
library(FitAR)
library(strucchange)
vcov.ri <- function(x, ...) kernHAC(x, kernel = "Quadratic Spectral",
prewhite = 1, approx = "AR(1)", ...)
library(readxl)
raw <- read_excel("NDRisk.xlsx", range ='main!A1:H51')
summary(raw)## Year Events Casualties Individual risk
## Min. :1970 Min. : 47.5 Min. : 6361 Min. : 1.062
## 1st Qu.:1982 1st Qu.:147.5 1st Qu.: 10502 1st Qu.: 1.970
## Median :1994 Median :266.7 Median : 18435 Median : 3.438
## Mean :1994 Mean :238.9 Mean : 52888 Mean : 10.142
## 3rd Qu.:2007 3rd Qu.:323.0 3rd Qu.: 46614 3rd Qu.: 8.869
## Max. :2019 Max. :443.2 Max. :384593 Max. :103.928
## Property risk Financial risk Real Losses Real GNI
## Min. : 5.294 Min. : 3.415 Min. : 7.71 Min. :1.994e+13
## 1st Qu.: 17.347 1st Qu.:10.513 1st Qu.: 36.21 1st Qu.:2.920e+13
## Median : 23.952 Median :15.278 Median : 79.66 Median :4.156e+13
## Mean : 34.807 Mean :18.867 Mean :110.20 Mean :4.607e+13
## 3rd Qu.: 43.228 3rd Qu.:22.760 3rd Qu.:155.96 3rd Qu.:6.256e+13
## Max. :115.829 Max. :56.870 Max. :474.95 Max. :8.443e+13DAT <- ts(raw, start=c(1970), end=c(2019),frequency=1)2 Simple trend analysis
Trend analsyis of the Natural Disaster Impact Indicators.
X = subset(DAT[,"Events"], start=11)
plot(X,main="Natural Disaster count")
trend = 1:length(X)
model = lm(X ~trend)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.917 -20.750 1.542 20.745 55.417
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 144.6455 8.7719 16.49 <2e-16 ***
## trend 6.7708 0.3728 18.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.22 on 38 degrees of freedom
## Multiple R-squared: 0.8967, Adjusted R-squared: 0.894
## F-statistic: 329.8 on 1 and 38 DF, p-value: < 2.2e-16X = log(DAT[,"Casualties"])
plot(X,ylab="log scale",main="Natural Disaster Casualties")
trend = 1:length(X)
model = lm(X ~trend)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4478 -0.7737 -0.3259 0.5880 2.6065
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.258703 0.316905 32.37 <2e-16 ***
## trend -0.005297 0.010816 -0.49 0.627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.104 on 48 degrees of freedom
## Multiple R-squared: 0.004972, Adjusted R-squared: -0.01576
## F-statistic: 0.2398 on 1 and 48 DF, p-value: 0.6265X = DAT[,"Real Losses"]
plot(X,ylab="2019$bn",main="Natural Disaster Real Losses")
trend = 1:length(X)
model = lm(X ~trend)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -99.16 -36.91 -10.54 20.54 283.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -16.0488 20.5995 -0.779 0.44
## trend 4.9508 0.7031 7.042 6.34e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 71.74 on 48 degrees of freedom
## Multiple R-squared: 0.5081, Adjusted R-squared: 0.4979
## F-statistic: 49.59 on 1 and 48 DF, p-value: 6.342e-09Trend analsyis of the Natural Disaster Risk Indicators.
X = log(DAT[,"Individual risk"])
plot(X,ylab="Log of Fatalities per million population",main="Natural Disaster Individual Risk")
trend = 1:length(X)
model = lm(X ~trend)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4412 -0.7583 -0.3139 0.5936 2.6475
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.01639 0.31668 6.367 6.87e-08 ***
## trend -0.02022 0.01081 -1.871 0.0675 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.103 on 48 degrees of freedom
## Multiple R-squared: 0.06796, Adjusted R-squared: 0.04854
## F-statistic: 3.5 on 1 and 48 DF, p-value: 0.06747X = subset( log(DAT[,"Property risk"]) , start=11)
plot(X,ylab="Log of risk",main="Natural Disaster Property Risk")
trend = 1:length(X)
model = lm(X ~trend)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.27183 -0.38819 -0.09984 0.40539 1.28947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.615760 0.193755 18.662 <2e-16 ***
## trend -0.008059 0.008236 -0.979 0.334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6013 on 38 degrees of freedom
## Multiple R-squared: 0.02458, Adjusted R-squared: -0.00109
## F-statistic: 0.9575 on 1 and 38 DF, p-value: 0.334X = subset( log(DAT[,"Financial risk"]) , start=11)
plot(X,ylab="Log of risk",main="Natural Disaster Financial Risk")
trend = 1:length(X)
model = lm(X ~trend)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72929 -0.37873 -0.03741 0.34790 1.09827
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.553680 0.154156 16.566 < 2e-16 ***
## trend 0.019041 0.006552 2.906 0.00608 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4784 on 38 degrees of freedom
## Multiple R-squared: 0.1818, Adjusted R-squared: 0.1603
## F-statistic: 8.444 on 1 and 38 DF, p-value: 0.0060773 Real GNI
We first study the yearly world Gross National Income (computed in constant 2010 US$). As usual for monetary amounts, we work with the logarithm of the original variable which displays two visual signals of non stationarity, namely trending up and slowly decaying autocorrelation function (ACF).
X = log(DAT[,"Real GNI"])
plot(X)
Acf(X)
The ARMA model selection based on the AIC criterion suggests 6 lags. We test for the presence of a unit root with the augmented Dickey-Fuller test as well as the Phillips-Perron test; both give large p-values, meaning we cannot reject the null hypothesis that there is one unit root.
SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 6 -399.1533 -314.3250
## 2 4 -399.0093 -315.7453
## 3 5 -397.7097 -315.3969
## 4 7 -397.2553 -312.7293
## 5 8 -397.0187 -311.3573adf.test(X, alternative = "explosive", k=6)##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -3.1281, Lag order = 6, p-value = 0.8782
## alternative hypothesis: explosivepp.test(X, alternative = "explosive")##
## Phillips-Perron Unit Root Test
##
## data: X
## Dickey-Fuller Z(alpha) = -18.326, Truncation lag parameter = 3, p-value
## = 0.9265
## alternative hypothesis: explosiveThe next step is to take first difference in the X series, which yields the GNI growth rate; we observe two signals for stationarity, mean-reversion and decaying ACF. The ARMA model selection based on the AIC criterion suggests 2 lags. Because the growth series is not trending we use the “stationary” option in both unit root tests. The test statistics are highly negative so that the p-values are very small; we may thus reject the null hypothesis that there is one unit root in the GNI growth rate. In other words, it is stationary meaning that GNI is integrated of order one (aka DS). We finally apply the ARMA model to the GNI growth; the AIC picks the AR(2) model (among meaningful values). The intercept value of 0.03 means that GNI is a random walk with drift growing at about 3% per year.
Y = diff(X)
plot(Y)
Acf(Y)
SelectModel(Y, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 2 -421.7692 -1.0617323
## 2 4 -420.1037 0.2175945
## 3 3 -419.8361 -1.5314104
## 4 1 -418.6879 -2.9920052
## 5 5 -418.3301 2.1469850adf.test(Y, alternative = "stationary" , k=2)## Warning in adf.test(Y, alternative = "stationary", k = 2): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Y
## Dickey-Fuller = -5.014, Lag order = 2, p-value = 0.01
## alternative hypothesis: stationarypp.test(Y, alternative = "stationary")## Warning in pp.test(Y, alternative = "stationary"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: Y
## Dickey-Fuller Z(alpha) = -30.288, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationaryarima(x = Y, order = c(2, 0, 0))##
## Call:
## arima(x = Y, order = c(2, 0, 0))
##
## Coefficients:
## ar1 ar2 intercept
## 0.3487 -0.3199 0.0293
## s.e. 0.1353 0.1374 0.0019
##
## sigma^2 estimated as 0.0001607: log likelihood = 144.36, aic = -280.72A last check is to detrend X to check whether it might be TS. It comes very close since the R-squared is exceptionally high.
trend = 1:length(X)
model = lm(X ~trend)
Z = ts(model$res,start=c(1970), end=c(2019),frequency=1)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.048090 -0.012142 -0.003926 0.009947 0.041683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.064e+01 5.458e-03 5614.8 <2e-16 ***
## trend 2.876e-02 1.863e-04 154.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01901 on 48 degrees of freedom
## Multiple R-squared: 0.998, Adjusted R-squared: 0.9979
## F-statistic: 2.385e+04 on 1 and 48 DF, p-value: < 2.2e-164 Disaster Real Losses
Real economic losses from natural disasters display two visual signals of non stationarity, namely trending up and slowly decaying autocorrelation function (ACF). The ARMA model selection based on the AIC criterion suggest 3 lags. We test for a unit root with Dickey-Fuller and Phillips-Perron; both give large p-values, meaning we cannot reject the null hypothesis that there is one unit root.
X = log(DAT[,"Real Losses"])
plot(X)
Acf(X)
SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 3 -45.19971 -46.12002
## 2 4 -44.49570 -44.34423
## 3 5 -42.71901 -44.15178
## 4 6 -42.38881 -43.27076
## 5 7 -41.43521 -41.78289adf.test(X, alternative = "explosive", k=3)##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -2.8313, Lag order = 3, p-value = 0.7592
## alternative hypothesis: explosivepp.test(X, alternative = "explosive")## Warning in pp.test(X, alternative = "explosive"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: X
## Dickey-Fuller Z(alpha) = -42.92, Truncation lag parameter = 3, p-value
## = 0.99
## alternative hypothesis: explosiveDue to the highly erratic behavior of real losses, the customary transformation is to detrend the series first. The fitted equation has highly significant parameters and features an 8% rate of growth (in real terms).
trend = 1:length(X)
model = lm(X ~trend)
Z = ts(model$res, start=c(1970), end=c(2019),frequency=1)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82962 -0.38553 -0.04132 0.41015 1.27233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.753377 0.152680 18.03 < 2e-16 ***
## trend 0.059404 0.005211 11.40 2.93e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5317 on 48 degrees of freedom
## Multiple R-squared: 0.7303, Adjusted R-squared: 0.7247
## F-statistic: 130 on 1 and 48 DF, p-value: 2.925e-15plot(Z)
Acf(Z)
The AIC suggests one lag or none at all. In both cases, the unit root tests return very low p-values; we may thus reject the null hypothesis that there is one unit root, meaning the X series is TS.
SelectModel(Z, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 0 -63.20317 3.463822
## 2 1 -61.72864 5.002933
## 3 2 -60.17151 5.392098
## 4 3 -59.68712 7.249289
## 5 4 -57.81844 9.248467adf.test(Z, alternative = "stationary", k = 0)## Warning in adf.test(Z, alternative = "stationary", k = 0): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Z
## Dickey-Fuller = -6.0087, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationaryadf.test(Z, alternative = "stationary", k = 1)## Warning in adf.test(Z, alternative = "stationary", k = 1): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Z
## Dickey-Fuller = -4.9062, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationarypp.test(Z, alternative = "stationary")## Warning in pp.test(Z, alternative = "stationary"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: Z
## Dickey-Fuller Z(alpha) = -42.92, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationaryWe lastly test for possible breaks in the trend. For that task, we must study the first difference and apply a test for structural change. The Bayesian Information Criteria recommeneded by Bai & Perron gives credence to an absence of break.
U = diff(X)
bp.ri <- breakpoints(U ~ 1)
summary(bp.ri) ##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = U ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 41
## m = 2 31 41
## m = 3 24 31 41
## m = 4 14 25 32 41
## m = 5 10 17 25 32 41
##
## Corresponding to breakdates:
##
## m = 1 2011
## m = 2 2001 2011
## m = 3 1994 2001 2011
## m = 4 1984 1995 2002 2011
## m = 5 1980 1987 1995 2002 2011
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 23.96 23.64 23.42 23.15 22.89 22.73
## BIC 111.78 118.91 126.24 133.45 140.68 148.11plot(bp.ri)
5 Individual Risk
Since the individual risk is also a ratio expressed in basis points, it is not transformed prior to the stationarity study; the temporal plot fails to display a clear trend and we further observe a rapidly decaying ACF. The ARMA model selection based on the AIC criterion suggest 6 lags. The Dickey-Fuller and Phillips-Perron unit root tests yield large p-values, leading to the conclusion that property risk is non-stationary.
X <- DAT[,"Individual risk"]
plot(X)
Acf(X)
SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 6 291.8304 0.8422301
## 2 8 293.5498 2.4141095
## 3 7 293.5922 0.5347569
## 4 0 293.7854 3.8937229
## 5 9 295.4136 3.5306689adf.test(X, alternative = "explosive", k=6)##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -2.8308, Lag order = 6, p-value = 0.759
## alternative hypothesis: explosiveadf.test(X, alternative = "explosive", k=0)## Warning in adf.test(X, alternative = "explosive", k = 0): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -10.065, Lag order = 0, p-value = 0.99
## alternative hypothesis: explosivepp.test(X, alternative = "explosive")## Warning in pp.test(X, alternative = "explosive"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: X
## Dickey-Fuller Z(alpha) = -52.074, Truncation lag parameter = 3, p-value
## = 0.99
## alternative hypothesis: explosiveThe search for a trend is unsucessful as the linear model returns a very low R-squared.
trend = 1:length(X)
model = lm(X ~trend)
Z = ts(model$res, start=c(1970), end=c(2019),frequency=1)
plot(Z)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.201 -9.576 -3.927 -1.794 86.891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.3185 5.2898 3.274 0.00197 **
## trend -0.2814 0.1805 -1.559 0.12560
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.42 on 48 degrees of freedom
## Multiple R-squared: 0.04819, Adjusted R-squared: 0.02836
## F-statistic: 2.43 on 1 and 48 DF, p-value: 0.1256The next transformation is the first difference. The ARMA model selection based on the AIC criterion suggest 5 lags. We test for a unit root with Dickey-Fuller and Phillips-Perron; both return very low p-values leading us to reject the null hypothesis that there is one unit root, meaning the individual risk series is DS. The suggested ARMA model has a significative intercept value of -0.57 which means that individual risk is expected to fall but importantly, as all autoregressive coefficients are negative, the level achieved in a given year is likely to be opposed to what happened in the 5 previous ones.
Y = diff(X)
plot(Y)
Acf(Y)
SelectModel(Y, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 5 293.1677 -12.653082
## 2 7 293.9119 -10.593158
## 3 6 294.3428 -12.592747
## 4 8 295.8554 -9.420050
## 5 9 296.7078 -8.420255adf.test(Y, alternative = "stationary" , k=5)## Warning in adf.test(Y, alternative = "stationary", k = 5): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Y
## Dickey-Fuller = -4.3738, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationarypp.test(Y, alternative = "stationary")## Warning in pp.test(Y, alternative = "stationary"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: Y
## Dickey-Fuller Z(alpha) = -58.901, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationaryarima(x = Y, order = c(5, 0, 0))##
## Call:
## arima(x = Y, order = c(5, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 intercept
## -0.9428 -0.7677 -0.7615 -0.7582 -0.6926 -0.5424
## s.e. 0.1480 0.1932 0.1868 0.1703 0.1246 0.4808
##
## sigma^2 estimated as 245.4: log likelihood = -206.44, aic = 426.87Lastly, we look at the possibility of structural breaks in the series. There might have been one in 1977, separating the deadly 1970s from the modern safer epoch.
U = X # or detrended Z or difference Y
bp.ri <- breakpoints(U ~ 1)
summary(bp.ri) ##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = U ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 7
## m = 2 7 41
## m = 3 7 33 41
## m = 4 7 22 33 41
## m = 5 7 15 22 33 41
## m = 6 7 14 21 28 35 42
##
## Corresponding to breakdates:
##
## m = 1 1976
## m = 2 1976 2010
## m = 3 1976 2002 2010
## m = 4 1976 1991 2002 2010
## m = 5 1976 1984 1991 2002 2010
## m = 6 1976 1983 1990 1997 2004 2011
##
## Fit:
##
## m 0 1 2 3 4 5 6
## RSS 17115.4 14217.7 13918.5 12789.3 12753.9 12626.2 13384.1
## BIC 441.5 440.1 446.8 450.4 458.1 465.4 476.2plot(bp.ri)
confid <- confint(bp.ri)
plot(U)
lines(bp.ri)
lines(confid)
coef(bp.ri, breaks = 1)## (Intercept)
## 1970 - 1976 29.009882
## 1977 - 2019 7.070416sapply(vcov(bp.ri, breaks = 1, vcov = vcov.ri), sqrt)## 1970 - 1976 1977 - 2019
## 10.405234 1.575265confint(bp.ri, breaks = 1, vcov = vcov.ri)##
## Confidence intervals for breakpoints
## of optimal 2-segment partition:
##
## Call:
## confint.breakpointsfull(object = bp.ri, breaks = 1, vcov. = vcov.ri)
##
## Breakpoints at observation number:
## 2.5 % breakpoints 97.5 %
## 1 5 7 26
##
## Corresponding to breakdates:
## 2.5 % breakpoints 97.5 %
## 1 1974 1976 1995fac.ri <- breakfactor(bp.ri, breaks = 1, label = "PRisk")
fm.ri <- lm(U ~ 0 + fac.ri)
plot(U)
lines(fitted(bp.ri, breaks = 1), col = 2)
summary(fm.ri)##
## Call:
## lm(formula = U ~ 0 + fac.ri)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.018 -5.567 -4.159 1.126 74.918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## fac.riPRisk1 29.010 6.505 4.460 4.94e-05 ***
## fac.riPRisk2 7.070 2.625 2.694 0.0097 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.21 on 48 degrees of freedom
## Multiple R-squared: 0.3612, Adjusted R-squared: 0.3346
## F-statistic: 13.57 on 2 and 48 DF, p-value: 2.129e-056 Property Risk
Since property risk is a ratio expressed in basis points, it is not transformed prior to the stationarity study; the temporal plot fails to display a clear trend and we further observe a rapidly decaying ACF. The ARMA model selection based on the AIC criterion suggest 2 lags. The two unit root tests previously employed yield large p-values, leading to the conclusion that property risk is non-stationary.
X <- DAT[,"Property risk"]
plot(X)
Acf(X)
SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 2 327.7857 0.4164036
## 2 3 328.5698 2.4152334
## 3 1 328.7379 -0.2908179
## 4 5 329.1622 2.1718711
## 5 0 329.8907 0.7832466adf.test(X, alternative = "explosive", k=2)##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -2.3674, Lag order = 2, p-value = 0.5731
## alternative hypothesis: explosivepp.test(X, alternative = "explosive")## Warning in pp.test(X, alternative = "explosive"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: X
## Dickey-Fuller Z(alpha) = -40.639, Truncation lag parameter = 3, p-value
## = 0.99
## alternative hypothesis: explosiveDue to the highly erratic behavior of this ratio, the customary transformation is to detrend the series first. The fitted equation has a very low R-squared and a non significant trend parameter, leading us to reject the TS hypothesis.
trend = 1:length(X)
model = lm(X ~trend)
Z = ts(model$res, start=c(1970), end=c(2019),frequency=1)
plot(Z)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.405 -17.316 -10.573 7.436 80.546
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.3392 7.7585 4.039 0.000192 ***
## trend 0.1360 0.2648 0.514 0.609945
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.02 on 48 degrees of freedom
## Multiple R-squared: 0.005464, Adjusted R-squared: -0.01526
## F-statistic: 0.2637 on 1 and 48 DF, p-value: 0.6099The logical next transformation is to apply “first difference.” The ARMA model selection based on the AIC criterion suggest 4 lags. We test for a unit root with Dickey-Fuller and Phillips-Perron; both return very low p-values leading us to reject the null hypothesis that there is a unit root. Hence, we conclude property risk to be DS. The suggested ARMA model has an intercept value of 0.09 which means that property risk is expected to vary very little but importantly, as all autoregressive coefficients are negative, the level achieved in a given year is likely to be opposed to what happened in the 4 previous ones; a phenomenon known as mean-reversion.
Y = diff(X)
plot(Y)
Acf(Y)
SelectModel(Y, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 4 324.0971 -16.80140
## 2 5 326.0395 -15.03352
## 3 2 326.6720 -14.42977
## 4 6 327.8346 -13.58215
## 5 3 327.9259 -18.74093adf.test(Y, alternative = "stationary" , k=4)## Warning in adf.test(Y, alternative = "stationary", k = 4): p-value smaller than
## printed p-value##
## Augmented Dickey-Fuller Test
##
## data: Y
## Dickey-Fuller = -4.2721, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationarypp.test(Y, alternative = "stationary")## Warning in pp.test(Y, alternative = "stationary"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: Y
## Dickey-Fuller Z(alpha) = -64.211, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationaryarima(x = Y, order = c(4, 0, 0))##
## Call:
## arima(x = Y, order = c(4, 0, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4 intercept
## -0.7721 -0.5324 -0.3645 -0.3233 0.0901
## s.e. 0.1329 0.1645 0.1616 0.1286 1.2001
##
## sigma^2 estimated as 595.4: log likelihood = -226.58, aic = 465.15Lastly, we look at the possibility of structural breaks in the series. The most likely ones are 1990 and 1999 with the in-between period of much higher risk (70 bp) when compared to the early and late periods (25 & 28 bp).
U = X # or detrended Z or difference Y
bp.ri <- breakpoints(U ~ 1)
summary(bp.ri) ##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = U ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 18
## m = 2 21 30
## m = 3 16 23 30
## m = 4 16 23 30 42
## m = 5 7 16 23 30 42
## m = 6 7 14 21 29 36 43
##
## Corresponding to breakdates:
##
## m = 1 1987
## m = 2 1990 1999
## m = 3 1985 1992 1999
## m = 4 1985 1992 1999 2011
## m = 5 1976 1985 1992 1999 2011
## m = 6 1976 1983 1990 1998 2005 2012
##
## Fit:
##
## m 0 1 2 3 4 5 6
## RSS 35236.5 30860.1 21132.3 19543.5 18868.1 18458.0 20072.7
## BIC 477.6 478.8 467.7 471.6 477.7 484.4 496.4plot(bp.ri)
confid <- confint(bp.ri)
plot(U)
lines(bp.ri)
lines(confid)
coef(bp.ri, breaks = 2)## (Intercept)
## 1970 - 1990 25.44150
## 1991 - 1999 70.52684
## 2000 - 2019 28.56589sapply(vcov(bp.ri, breaks = 2, vcov = vcov.ri), sqrt)## 1970 - 1990 1991 - 1999 2000 - 2019
## 4.586194 5.535064 2.954736confint(bp.ri, breaks = 2, vcov = vcov.ri)##
## Confidence intervals for breakpoints
## of optimal 3-segment partition:
##
## Call:
## confint.breakpointsfull(object = bp.ri, breaks = 2, vcov. = vcov.ri)
##
## Breakpoints at observation number:
## 2.5 % breakpoints 97.5 %
## 1 19 21 24
## 2 28 30 32
##
## Corresponding to breakdates:
## 2.5 % breakpoints 97.5 %
## 1 1988 1990 1993
## 2 1997 1999 2001fac.ri <- breakfactor(bp.ri, breaks = 2, label = "PRisk")
fm.ri <- lm(U ~ 0 + fac.ri)
plot(U)
lines(fitted(bp.ri, breaks = 2), col = 2)
summary(fm.ri)##
## Call:
## lm(formula = U ~ 0 + fac.ri)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.577 -11.212 -5.575 7.122 78.398
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## fac.riPRisk1 25.442 4.627 5.498 1.53e-06 ***
## fac.riPRisk2 70.527 7.068 9.978 3.44e-13 ***
## fac.riPRisk3 28.566 4.741 6.025 2.47e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.2 on 47 degrees of freedom
## Multiple R-squared: 0.7794, Adjusted R-squared: 0.7654
## F-statistic: 55.36 on 3 and 47 DF, p-value: 1.846e-157 Financial Risk
Since financial risk is a ratio expressed in basis points, it is not transformed prior to the stationarity study; the temporal plot displays a positive trend and we further observe a slowly decaying ACF. The ARMA model selection based on the AIC criterion suggest 6 lags. The two unit root tests yield large p-values, leading to the conclusion that financial risk is non-stationary.
X <- DAT[,"Financial risk"]
plot(X)
Acf(X)
SelectModel(X, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 6 244.4509 -7.912183
## 2 7 246.4333 -5.952399
## 3 3 248.1233 -3.152902
## 4 8 248.3983 -4.040886
## 5 5 248.8851 -9.889972adf.test(X, alternative = "explosive", k=6)##
## Augmented Dickey-Fuller Test
##
## data: X
## Dickey-Fuller = -1.8177, Lag order = 6, p-value = 0.3526
## alternative hypothesis: explosivepp.test(X, alternative = "explosive")## Warning in pp.test(X, alternative = "explosive"): p-value smaller than printed
## p-value##
## Phillips-Perron Unit Root Test
##
## data: X
## Dickey-Fuller Z(alpha) = -40.305, Truncation lag parameter = 3, p-value
## = 0.99
## alternative hypothesis: explosiveDue to the highly erratic behavior of this ratio, the customary transformation is to detrend the series first. The fitted model has a highly significant positive slope parameter of nearly half a basis point per year.
trend = 1:length(X)
model = lm(X ~trend)
Z = ts(model$res, start=c(1970), end=c(2019),frequency=1)
summary(model)##
## Call:
## lm(formula = X ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.429 -6.101 -2.747 4.603 33.173
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.6708 3.0588 2.181 0.0341 *
## trend 0.4783 0.1044 4.581 3.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.65 on 48 degrees of freedom
## Multiple R-squared: 0.3042, Adjusted R-squared: 0.2897
## F-statistic: 20.99 on 1 and 48 DF, p-value: 3.299e-05plot(Z)
Acf(Z)
The AIC suggests 2 lags. The unit root tests return very low p-values; we may thus reject the null hypothesis that there is one unit root, meaning that the financial risk series is TS.
SelectModel(Z, ARModel="AR", Criterion="AIC", Candidates = 15, Best=5)## p AIC-Exact AIC-Approx
## 1 2 234.9783 1.4333137
## 2 4 235.3328 2.1488828
## 3 3 236.2428 0.2803927
## 4 0 236.5397 3.9000712
## 5 6 236.6140 2.9748952adf.test(Z, alternative = c("stationary"), k = 2)##
## Augmented Dickey-Fuller Test
##
## data: Z
## Dickey-Fuller = -3.9464, Lag order = 2, p-value = 0.01904
## alternative hypothesis: stationarypp.test(Z, alternative = c("stationary"))## Warning in pp.test(Z, alternative = c("stationary")): p-value smaller than
## printed p-value##
## Phillips-Perron Unit Root Test
##
## data: Z
## Dickey-Fuller Z(alpha) = -40.305, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationaryWe lastly test for possible breaks in the trend and study the first difference. We apply the previous test for structural change. The Bayesian Information Criteria gives credence to an absence of break.
U = diff(X)
bp.ri <- breakpoints(U ~ 1)
plot(bp.ri)
summary(bp.ri) ##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = U ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 41
## m = 2 31 41
## m = 3 25 32 41
## m = 4 17 25 32 41
## m = 5 10 17 25 32 41
##
## Corresponding to breakdates:
##
## m = 1 2011
## m = 2 2001 2011
## m = 3 1995 2002 2011
## m = 4 1987 1995 2002 2011
## m = 5 1980 1987 1995 2002 2011
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 10266.8 10026.7 9880.6 9690.2 9557.7 9521.1
## BIC 408.7 415.4 422.4 429.3 436.4 444.0End of File