Timberland Project Report
Yijia Yu
18/06/2020
1. Introduction
Timberland investment instruments are often used as an alternative asset class in the public and private pension fund maneger’s portfolio. Unlike the traditional financial tools such as stocks, bonds, cash, timberland asset has the benefit of low correlation to the stock market and the tendency of rise over in a long run. Therefore, timberland assets not only makes the portfolio more diversified but also minimizes the volatility. A powerful component of timberland return drivers is biological growth, which is independent of all factors typically impacting other investments and makes timberland asset distinguishing from others.
Timberland Industry in Canada
The forestry industry is a key contributor to Canada’s economy, with real GDP growth recorded at 0.7 percent in 2017. This comes as no surprise, as there were around 308.94 million hectares of forest area recorded in Canada in 2019
The majority of Canada’s forests are owned by the provinces in which they are located, with just 1.6 percent being owned federally. The provincial and territorial governments are responsible for managing forest conservation, with strict laws in place to shape forest practices. These stringent laws are particularly important when it comes to halting logging practices and encouraging reforestation, as extensive deforestation was carried out by European settlers in Canada throughout the 18th and 19th centuries. The wood product manufacturing sector is the most lucrative of the forestry industries, with a nominal GDP of approximately 11 billion Canadian dollars in 2017. Western Forest Products, a lumber company based in British Columbia, is one of the key players in the forestry industry with an annual revenue of around 1.2 billion Canadian dollars in 2018. As a result, British Columbia employs by far the most people in the forestry and logging industry, followed by Quebec. The United States is the main trading partner of Canada when it comes to forest resources, with Canadian exports worth around 24.24 billion Canadian dollars in 2017. Lumber, sawmill and millwork products are the most traded Canadian forest products, with exports worth just over 15 billion Canadian dollars in 2018 globally.
In spite of today’s challenging economic environment, the long-term supply and demand fundamentals of timberland bode well for the patient investor who thinks in terms of decades and not quarters. A growing canadian and global population with a rising standard of living will continue to increase demand on a finite supply of timberland.
About this study
The main purpose of this study is to investigate long-term financial performance metrics of both private and public timberland investments. I will employ the modern portfolio theory to evaluate its performance and to analyze its diversification effects.
Data used in this study
Private timberland market Returns of private equity timberland investment performance is approximated by the NCREIF Timberland Index (NTI), which is a quarterly time series composite return measure of investment performance of a large pool of individual U.S. timber properties acquired in the private market for investment purposes only.
Public timberland market As of July 2019, there is four publicly traded REITs that specialize in timberland. Together, they own more than 17 million acres of timberland properties in the U.S., Canada, New Zealand and other locations. The four REITs compose the timber sector of the FTSE NAREIT All REIT index are listed below:
- Potlatch
- CatchMark
- Rayonier
- Weyerhaeuser
2. Methods
2.1 Build Portfoilos
To examine the role of timberland assets in a mixed portfolio, assets including large-cap stocks, small-cap stocks, treasury bonds, and treasury bills are considered in this study. Returns on those asset classes are approximated by the S&P 500 Index (SP500), Russell 2000 Index (RU2000), 10-Year Treasury bonds yield, and 3-month treasury bills.
2.2 Asset return volatility
Generalized autoregressive conditional heteroskedasticity (GARCH) models are used to exam the time-varying conditional variances.
Results
3.1 Descriptive analysis
This is a interactively graphical representation of selected quarterly return indices from 1987 Q1 to 2020 Q1. Some clear pattern could be observed. First, NCREIF is much less volatile than REIT, SP500, RU2000. Second, REIT follows SP500 closely but NTI doesn’t. For example, both REIT and SP500 suffered from a big loss during 2008 financial crisis, whereas NCREIF kept a positive returns until 2009/Q2, where they only got a slightly negative return value. Some people said it’s because timberland market tends to lag the public equity market. But some also argues that it’s because of the biological growth of trees wouldn’t stop as the market went into depression.
dygraph(allxts*100, main = "Quarterly returns of analyzed assets", ylab = "Return%")%>%
dyHighlight(highlightSeriesOpts = list(strokeWidth = 3))%>%
dyRangeSelector()%>%
dyLegend(width = 500)%>%
dyOptions(colors = RColorBrewer::brewer.pal(7, "Dark2"))The unconditional correlations of the return series were studied by calculating the correlation coefficients and by plotting the pairwise scatter plots, as showed in the following graph. In this portfolio, the correlation of NCREIF with all financial assets are uniformly low, except for the 3 month Treasury bill (CA). The correlation between big cap (SP500) and small cap (RU2000) are highly significant, and they’re almost linearly correlated. Unlike private timber market, the timber REITs shows a closely correlation with the stock market and 10 year treasury bonds. Returns on most asset classes are positively skewed. However, by using the Jarque-Bera test statistic, most of asset returns are not normally distributed.
ggpairs(as.data.frame(allts))+theme_bw()
# normality
apply(allts,2,jarque.bera.test)## $NTI
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 365.15, df = 2, p-value < 2.2e-16
##
##
## $NPI
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 415.78, df = 2, p-value < 2.2e-16
##
##
## $SP500
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 20.6, df = 2, p-value = 3.364e-05
##
##
## $RU2000
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 17.542, df = 2, p-value = 0.0001551
##
##
## $`T-bonds10Y`
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 1.3328, df = 2, p-value = 0.5136
##
##
## $`T-bills3M`
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 27.516, df = 2, p-value = 1.059e-06
##
##
## $`Timber REITs`
##
## Jarque Bera Test
##
## data: newX[, i]
## X-squared = 2400.2, df = 2, p-value < 2.2e-16Summary statistics of all asset classes from 1987 Q1 to 2010 Q1 are reported in the following table. On a quarterly basis, NTI (private timberland index) has the highest mean return of 2.71% and a moderate standard deviation of 3.78%. Timber REITs have a highest standard deviation but a relative lower mean of 2.06%. Treasury bills have both lowest return and standard deviation.
describe(allxts*100)[,-1]## n mean sd median trimmed mad min max range skew
## NTI 133 2.71 3.78 1.53 2.07 1.51 -6.54 22.34 28.88 2.13
## NPI 133 1.90 2.09 2.11 2.18 1.13 -8.29 5.43 13.72 -2.23
## SP500 133 2.13 8.02 2.93 2.72 5.22 -23.23 20.87 44.09 -0.76
## RU2000 133 2.21 10.60 3.26 2.90 8.18 -30.89 29.16 60.05 -0.65
## T-bonds10Y 133 1.54 2.30 1.48 1.51 2.50 -3.39 8.04 11.42 0.14
## T-bills3M 133 0.04 0.03 0.03 0.03 0.03 0.00 0.14 0.13 1.07
## Timber REITs 133 2.06 11.46 2.55 2.04 7.42 -33.14 85.45 118.58 2.47
## kurtosis se
## NTI 6.73 0.33
## NPI 7.24 0.18
## SP500 1.11 0.70
## RU2000 1.13 0.92
## T-bonds10Y -0.43 0.20
## T-bills3M 0.44 0.00
## Timber REITs 19.85 0.99Three risk measures are applied to individual asset, and the results are reported in the following table.
res = rbind(
apply(allts, 2, sd),
VaR(allts, p=.95, method="modified"),
ES(allts, p=.95, method="modified")
)
row.names(res)[1] = "SD"
row.names(res)[3] = "CVaR"
res## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## SD 0.037831038 0.02091420 0.08019994 0.1059724 0.02296897 3.349197e-04
## VaR -0.003228681 -0.02349186 -0.12481671 -0.1679863 -0.02142668 -4.166894e-05
## CVaR -0.107178392 -0.05549443 -0.18308384 -0.2478092 -0.02867966 -2.001534e-04
## Timber REITs
## SD 0.11455577
## VaR -0.02612674
## CVaR -0.18075306# stationarity
apply(allts,2,adf.test)## $NTI
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -3.9265, Lag order = 5, p-value = 0.01501
## alternative hypothesis: stationary
##
##
## $NPI
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -4.2311, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $SP500
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -4.0533, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $RU2000
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -5.1613, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $`T-bonds10Y`
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -5.3722, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $`T-bills3M`
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -3.1882, Lag order = 5, p-value = 0.09292
## alternative hypothesis: stationary
##
##
## $`Timber REITs`
##
## Augmented Dickey-Fuller Test
##
## data: newX[, i]
## Dickey-Fuller = -4.9764, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryapply(allts,2,pp.test)## $NTI
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -177.73, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
##
##
## $NPI
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -29.936, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
##
##
## $SP500
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -135.84, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
##
##
## $RU2000
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -137.45, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
##
##
## $`T-bonds10Y`
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -122.7, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary
##
##
## $`T-bills3M`
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -16.293, Truncation lag parameter = 4, p-value
## = 0.1758
## alternative hypothesis: stationary
##
##
## $`Timber REITs`
##
## Phillips-Perron Unit Root Test
##
## data: newX[, i]
## Dickey-Fuller Z(alpha) = -134.5, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryop = par(mfrow=c(1,2))
#test stationarity and normality
for(i in 1:ncol(allxts)){
#plot acf
acf(na.omit(allxts[,i])^2,lag.max = 100,xlab = colnames(allxts)[i], ylab = 'ACF', main=' ')
#plot pacf
pacf(na.omit(allxts[,i])^2,lag.max = 100,xlab = colnames(allxts)[i], ylab = 'PACF', main=' ')
}
3.2 Portfolio Optimization
library(RColorBrewer)
library(fPortfolio)## Loading required package: timeDate##
## Attaching package: 'timeDate'## The following objects are masked from 'package:PerformanceAnalytics':
##
## kurtosis, skewness## Loading required package: timeSeries##
## Attaching package: 'timeSeries'## The following object is masked from 'package:psych':
##
## outlier## The following object is masked from 'package:zoo':
##
## time<-## Loading required package: fBasics##
## Attaching package: 'fBasics'## The following object is masked from 'package:psych':
##
## tr## Loading required package: fAssetspacman::p_load(matrixcalc,knitr,dygraphs,ggthemes,highcharter,viridis,tibbletime,timetk,tidyquant,tidyverse,fPortfolio,xts)
#0. Prepare data
returns.data <- na.omit(allts)
returns.data %>%tk_tbl()%>%head()## # A tibble: 6 x 8
## index NTI NPI SP500 RU2000 `T-bonds10Y` `T-bills3M` `Timber REITs`
## <year> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1987 … -6.00e-4 0.0183 0.205 0.238 0.0148 0.000729 0.252
## 2 1987 … 4.91e-2 0.0119 0.0422 -0.0113 -0.0189 0.000813 -0.0558
## 3 1987 … 1.19e-1 0.0209 0.0587 0.0383 -0.0292 0.000905 -0.00576
## 4 1987 … 7.85e-2 0.0267 -0.232 -0.294 0.0584 0.000835 -0.0944
## 5 1988 … 3.98e-2 0.0184 0.0478 0.184 0.0358 0.000832 0.0147
## 6 1988 … 6.33e-2 0.02 0.0564 0.0615 0.00986 0.000907 0.0210returns_ts <- as.timeSeries(returns.data)
spec <- portfolioSpec()
setSolver(spec) <- "solveRquadprog"
setNFrontierPoints(spec) <-25
#1. calculating equal weight portfolio
nAssets <- ncol(returns_ts)
setWeights(spec)<-rep(1/nAssets, times = nAssets)
constraints <- 'LongOnly'
portfolioConstraints(returns_ts, spec, constraints)##
## Title:
## Portfolio Constraints
##
## Lower/Upper Bounds:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M Timber REITs
## Lower 0 0 0 0 0 0 0
## Upper 1 1 1 1 1 1 1
##
## Equal Matrix Constraints:
## ceq NTI NPI SP500 RU2000 T-bonds10Y T-bills3M Timber REITs
## Budget -1 -1 -1 -1 -1 -1 -1 -1
## attr(,"na.action")
## Return
## 1
## attr(,"class")
## [1] "omit"
##
## Cardinality Constraints:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M Timber REITs
## Lower 0 0 0 0 0 0 0
## Upper 1 1 1 1 1 1 1# calculate the properties of the portfolio
# Now let us display the results from the equal weights portfolio, the assignment of weights, and the attribution of returns and risk.
ewPortfolio <- feasiblePortfolio( returns_ts, spec, constraints)
print(ewPortfolio)##
## Title:
## MV Feasible Portfolio
## Estimator: covEstimator
## Solver: solveRquadprog
## Optimize: minRisk
## Constraints: LongOnly
##
## Portfolio Weights:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## 0.1429 0.1429 0.1429 0.1429 0.1429 0.1429
## Timber REITs
## 0.1429
##
## Covariance Risk Budgets:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## 0.0175 0.0124 0.2697 0.3675 -0.0120 0.0002
## Timber REITs
## 0.3447
##
## Target Returns and Risks:
## mean Cov CVaR VaR
## 0.0180 0.0375 0.0807 0.0624
##
## Description:
## Thu Jun 18 15:41:49 2020 by user:col <- divPalette(ncol(returns_ts), "RdBu")
weightsPie(ewPortfolio, radius = 0.7, col = col)
mtext(text = "Equally Weighted MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
weightedReturnsPie(ewPortfolio, radius = 0.7, col = col)
mtext(text = "Equally Weighted MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
covRiskBudgetsPie(ewPortfolio, radius = 0.7, col = col)
mtext(text = "Equally Weighted MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
# COMPUTE A MINIMUM RISK EFFICIENT PORTFOLIO!
minriskSpec <- portfolioSpec()
targetReturn <- getTargetReturn(ewPortfolio@portfolio)["mean"]
setTargetReturn(minriskSpec) <- targetReturn
minriskPortfolio <- efficientPortfolio(
data = returns_ts,
spec = minriskSpec,
constraints = "LongOnly")
print(minriskPortfolio)##
## Title:
## MV Efficient Portfolio
## Estimator: covEstimator
## Solver: solveRquadprog
## Optimize: minRisk
## Constraints: LongOnly
##
## Portfolio Weights:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## 0.1551 0.4308 0.0000 0.0147 0.3134 0.0666
## Timber REITs
## 0.0194
##
## Covariance Risk Budgets:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## 0.2350 0.4561 0.0000 0.0181 0.2680 0.0005
## Timber REITs
## 0.0222
##
## Target Returns and Risks:
## mean Cov CVaR VaR
## 0.0180 0.0126 0.0141 0.0014
##
## Description:
## Thu Jun 18 15:41:49 2020 by user:col <- qualiPalette(ncol(returns_ts), "Dark2")
weightsPie(minriskPortfolio, radius = 0.7, col = col)
mtext(text = "Minimal Risk MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
weightedReturnsPie(minriskPortfolio, radius = 0.7, col = col)
mtext(text = "Minimal Risk MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
covRiskBudgetsPie(minriskPortfolio, radius = 0.7, col = col)
mtext(text = "Minimal Risk MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
# Tangency Portfolio
tgSpec <- portfolioSpec()
setRiskFreeRate(tgSpec) <- 0
tgPortfolio <- tangencyPortfolio(
data = returns_ts,
spec = tgSpec,
constraints = "LongOnly")
print(tgPortfolio)##
## Title:
## MV Tangency Portfolio
## Estimator: covEstimator
## Solver: solveRquadprog
## Optimize: minRisk
## Constraints: LongOnly
##
## Portfolio Weights:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## 0.0017 0.0161 0.0000 0.0008 0.0090 0.9725
## Timber REITs
## 0.0000
##
## Covariance Risk Budgets:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M
## 0.0504 0.3439 0.0000 0.0197 0.1553 0.4301
## Timber REITs
## 0.0006
##
## Target Returns and Risks:
## mean Cov CVaR VaR
## 9e-04 5e-04 5e-04 -1e-04
##
## Description:
## Thu Jun 18 15:41:50 2020 by user:col <- seqPalette(ncol(returns_ts), "BuPu")
weightsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
weightedReturnsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
covRiskBudgetsPie(tgPortfolio, box = FALSE, col = col)
mtext(text = "Tangency MV Portfolio", side = 3, line = 1.5,
font = 2, cex = 0.7, adj = 0)
#Now we are ready to start the optmization and pass data, specification and constraints to the function portfolioFrontier(). We can examine the result using a standard print() command.
# compute the efficient frontier
setNFrontierPoints(spec) <-25
frontier <- portfolioFrontier(returns_ts, spec, constraints)
print(frontier)##
## Title:
## MV Portfolio Frontier
## Estimator: covEstimator
## Solver: solveRquadprog
## Optimize: minRisk
## Constraints: LongOnly
## Portfolio Points: 5 of 25
##
## Portfolio Weights:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M Timber REITs
## 1 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000
## 7 0.0572 0.1661 0.0000 0.0058 0.1191 0.6448 0.0070
## 13 0.1172 0.3281 0.0000 0.0113 0.2381 0.2908 0.0146
## 19 0.2655 0.4697 0.0000 0.0167 0.2281 0.0000 0.0201
## 25 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
##
## Covariance Risk Budgets:
## NTI NPI SP500 RU2000 T-bonds10Y T-bills3M Timber REITs
## 1 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000
## 7 0.2248 0.4566 0.0000 0.0187 0.2645 0.0146 0.0209
## 13 0.2331 0.4563 0.0000 0.0182 0.2674 0.0031 0.0219
## 19 0.4607 0.3970 0.0000 0.0198 0.1021 0.0000 0.0204
## 25 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
##
## Target Returns and Risks:
## mean Cov CVaR VaR
## 1 0.0004 0.0003 0.0000 0.0000
## 7 0.0071 0.0048 0.0054 0.0004
## 13 0.0138 0.0096 0.0107 0.0010
## 19 0.0205 0.0149 0.0173 0.0005
## 25 0.0271 0.0378 0.0211 0.0026
##
## Description:
## Thu Jun 18 15:41:50 2020 by user:# plot efficient frontier
tailoredFrontierPlot(object = frontier, mText = "MV Portfolio - LongOnly Constraints",
risk = "Cov")
# plot weights
weightsPlot(frontier, col=rainbow(dim(returns.data)[2]))
text <- "Mean-Variance Portfolio - Long Only Constraints"
mtext(text, side = 3, line = 3, font = 2, cex = 0.9)
weightedReturnsPlot(frontier)
covRiskBudgetsPlot(frontier)
set.seed(12)
frontierPlot(object = frontier, pch = 19, cex = 0.5)
monteCarloPoints(object = frontier, mcSteps = 5000, pch = 19,
cex = 1,col="#D53E4F")
frontierpts <- frontierPoints(object = frontier)
lines(frontierpts, col = "blue", lwd = 2)
GARCH model estimation, Backtesting the risk model and Forecasting
#Ncreif
xts.nti <- allxts[, "NTI"]
xts.npi <- allxts[, "NPI"]
REITs <- allts[, "Timber REITs"]
# ARIMA
# ncreif
model.arima = auto.arima(xts.nti , trace = F , ic = 'bic')
model.arima## Series: xts.nti
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## -0.9517 0.0345 -0.6404
## s.e. 0.0473 0.1038 0.1170
##
## sigma^2 estimated as 0.001036: log likelihood=267.04
## AIC=-526.08 AICc=-525.77 BIC=-514.55checkresiduals(model.arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2)
## Q* = 7.0216, df = 5, p-value = 0.219
##
## Model df: 3. Total lags used: 8model.arima3 = auto.arima(xts.npi, stationary = T, trace = F , ic = 'bic')
model.arima3## Series: xts.npi
## ARIMA(2,0,2)(0,0,2)[4] with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 ma2 sma1 sma2 mean
## 1.4429 -0.7121 -0.6924 0.4407 0.5504 0.1771 0.0187
## s.e. 0.1290 0.1038 0.1408 0.0985 0.1120 0.0922 0.0044
##
## sigma^2 estimated as 0.0001201: log likelihood=413.99
## AIC=-811.97 AICc=-810.81 BIC=-788.85checkresiduals(model.arima3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,0,2)(0,0,2)[4] with non-zero mean
## Q* = 0.38528, df = 3, p-value = 0.9433
##
## Model df: 7. Total lags used: 10# REIT
model.arima2 = auto.arima(REITs, trace = T , ic = 'bic')##
## ARIMA(2,0,2)(1,0,1)[4] with non-zero mean : -164.7779
## ARIMA(0,0,0) with non-zero mean : -190.1259
## ARIMA(1,0,0)(1,0,0)[4] with non-zero mean : -181.7281
## ARIMA(0,0,1)(0,0,1)[4] with non-zero mean : -182.1282
## ARIMA(0,0,0) with zero mean : -190.7472
## ARIMA(0,0,0)(1,0,0)[4] with non-zero mean : -185.2885
## ARIMA(0,0,0)(0,0,1)[4] with non-zero mean : -185.2829
## ARIMA(0,0,0)(1,0,1)[4] with non-zero mean : -181.1828
## ARIMA(1,0,0) with non-zero mean : -186.6041
## ARIMA(0,0,1) with non-zero mean : -187.0044
## ARIMA(1,0,1) with non-zero mean : -183.2944
##
## Best model: ARIMA(0,0,0) with zero meancheckresiduals(model.arima2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,0) with zero mean
## Q* = 5.1745, df = 8, p-value = 0.7388
##
## Model df: 0. Total lags used: 8# box test
Box.test(model.arima$residuals^2,lag = 12,type="Ljung-Box")##
## Box-Ljung test
##
## data: model.arima$residuals^2
## X-squared = 24.212, df = 12, p-value = 0.01903Box.test(model.arima2$residuals^2, lag = 12, type="Ljung-Box")##
## Box-Ljung test
##
## data: model.arima2$residuals^2
## X-squared = 1.3134, df = 12, p-value = 0.9999Box.test(model.arima3$residuals^2, lag = 12, type="Ljung-Box")##
## Box-Ljung test
##
## data: model.arima3$residuals^2
## X-squared = 1.8949, df = 12, p-value = 0.9996#colnames(allts)
for(i in 1:ncol(allts)){
y = allts[,i]
res=Box.test((y - mean(y))^2, lag = 12,type="Ljung-Box")
print(res$p.value)
}## [1] 0.002018172
## [1] 1.529996e-10
## [1] 0.2494003
## [1] 0.855529
## [1] 0.8463122
## [1] 0
## [1] 0.9999845# GARCH model
library("rugarch")## Loading required package: parallel##
## Attaching package: 'rugarch'## The following object is masked from 'package:purrr':
##
## reduce## The following object is masked from 'package:stats':
##
## sigmaspec = ugarchspec(variance.model = list(model="sGARCH",garchOrder=c(1,1)),
mean.model = list(armaOrder=c(1,1)),
distribution.model = "std")
# Nelson's egarch model
egarch.spec = ugarchspec(variance.model=list(model="eGARCH",garchOrder=c(1,1)),
mean.model=list(armaOrder=c(0,0)),
distribution.model="std")
garch.fit1 = ugarchfit(egarch.spec, xts.nti)
# simulation
set.seed(123)
tseq = seq(as.Date("1987/1/1"), as.Date("2020/6/1"), "weeks")
sim = ugarchsim(garch.fit1,n.sim=length(tseq), n.start=0, m.sim=1, startMethod="sample")
simseries = xts(sim@simulation$seriesSim, order.by = tseq)
auto.arima(simseries , trace = T , ic = 'bic')##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,0,2) with non-zero mean : -11942.31
## ARIMA(0,0,0) with non-zero mean : -11971.15
## ARIMA(1,0,0) with non-zero mean : -11963.77
## ARIMA(0,0,1) with non-zero mean : -11964.52
## ARIMA(0,0,0) with zero mean : -10316.91
## ARIMA(1,0,1) with non-zero mean : -11956.67
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(0,0,0) with non-zero mean : -11971.15
##
## Best model: ARIMA(0,0,0) with non-zero mean## Series: simseries
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 0.0098
## s.e. 0.0002
##
## sigma^2 estimated as 6.067e-05: log likelihood=5993.04
## AIC=-11982.08 AICc=-11982.07 BIC=-11971.15garch.fit1 = ugarchfit(egarch.spec, xts.nti)
# simulation
set.seed(1234)
sim = ugarchsim(garch.fit1,n.sim=2000, n.start=0, m.sim=1, startMethod="sample")
simseries = as.data.frame(sim@simulation$seriesSim)
auto.arima(simseries, stationary = TRUE , trace = T , ic = 'bic')##
## Fitting models using approximations to speed things up...
##
## ARIMA(2,0,2) with non-zero mean : -14451.07
## ARIMA(0,0,0) with non-zero mean : -14480.05
## ARIMA(1,0,0) with non-zero mean : -14472.45
## ARIMA(0,0,1) with non-zero mean : -14472.7
## ARIMA(0,0,0) with zero mean : -12096.21
## ARIMA(1,0,1) with non-zero mean : -14464.86
##
## Now re-fitting the best model(s) without approximations...
##
## ARIMA(0,0,0) with non-zero mean : -14480.05
##
## Best model: ARIMA(0,0,0) with non-zero mean## Series: simseries
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 0.0098
## s.e. 0.0001
##
## sigma^2 estimated as 4.17e-05: log likelihood=7247.62
## AIC=-14491.25 AICc=-14491.24 BIC=-14480.05plot(garch.fit1, which="all")##
## please wait...calculating quantiles...
# backtesting model
garchroll1 <- ugarchroll(egarch.spec, data=simseries, n.start = 1000, refit.every = 100, refit.window = "moving",VaR.alpha = 0.01,solver="hybrid", fit.control = list(scale = 1))
report(garchroll1, type = "VaR",VaR.alpha = 0.01, conf.level = 0.99)## VaR Backtest Report
## ===========================================
## Model: eGARCH-std
## Backtest Length: 1000
## Data:
##
## ==========================================
## alpha: 1%
## Expected Exceed: 10
## Actual VaR Exceed: 13
## Actual %: 1.3%
##
## Unconditional Coverage (Kupiec)
## Null-Hypothesis: Correct Exceedances
## LR.uc Statistic: 0.831
## LR.uc Critical: 6.635
## LR.uc p-value: 0.362
## Reject Null: NO
##
## Conditional Coverage (Christoffersen)
## Null-Hypothesis: Correct Exceedances and
## Independence of Failures
## LR.cc Statistic: 2.833
## LR.cc Critical: 9.21
## LR.cc p-value: 0.243
## Reject Null: NOplot(garchroll1, which="all")
plot(garchroll1, which=4)
# Forecasting Risk and VaR
garchfcst <- ugarchforecast(garch.fit1, n.ahead = 12)
garchfcst##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: eGARCH
## Horizon: 12
## Roll Steps: 0
## Out of Sample: 0
##
## 0-roll forecast [T0=2020 Q1]:
## Series Sigma
## T+1 0.01004 0.010373
## T+2 0.01004 0.010281
## T+3 0.01004 0.010192
## T+4 0.01004 0.010107
## T+5 0.01004 0.010025
## T+6 0.01004 0.009946
## T+7 0.01004 0.009871
## T+8 0.01004 0.009799
## T+9 0.01004 0.009729
## T+10 0.01004 0.009662
## T+11 0.01004 0.009598
## T+12 0.01004 0.009536plot(garchfcst,which=1)
plot(garchfcst,which=3)
garch.fit2=ugarchfit(egarch.spec,data=xts.nti, solver="hybrid", out.sample=5)
garchfcst2<-ugarchforecast(garch.fit2, data = NULL, n.ahead = 10, n.roll = 5, external.forecasts = list(mregfor = NULL, vregfor = NULL))
garchfcst2##
## *------------------------------------*
## * GARCH Model Forecast *
## *------------------------------------*
## Model: eGARCH
## Horizon: 10
## Roll Steps: 5
## Out of Sample: 5
##
## 0-roll forecast [T0=2018 Q4]:
## Series Sigma
## T+1 0.01072 0.008570
## T+2 0.01072 0.008220
## T+3 0.01072 0.007884
## T+4 0.01072 0.007562
## T+5 0.01072 0.007254
## T+6 0.01072 0.006958
## T+7 0.01072 0.006675
## T+8 0.01072 0.006404
## T+9 0.01072 0.006143
## T+10 0.01072 0.005894plot(garchfcst2,which=2)
plot(garchfcst2,which=4)