计量经济学方法（An Econometric Approach）下的VaR计算
目录
1 基于innovation服从高斯分布下的算法——将ARMA序列转化为MA(∞)序列
1.1 unit root检验与处理
对于非平稳时间序列，我们可以从它的ACF图就能判断，即会存在很大的相关性。但判断其是否需要差分等处理则需要用更正式的统计检验。之所以要求序列平稳是为了防止出现伪估计等状况。
1.1.1 检验
此处方法较多，我们这里可采用ADF test来检验序列是否存在单位根和是否平稳。 若接受原假设，则存在单位根。若拒绝原假设，则序列平稳。
1.1.2 处理
若序列存在单位根，则进行差分处理，然后再进行单位根检验。
1.2平稳序列的相关性检验与建模（ARMA）
通过对处理后的平稳序列进行相关性检验，若存在ARMA，则建立相应 model，
1.2.1 相关性检验以及均值检验
1.图形检验
看了几篇文章，发现ACF在检验上为MA检验服务的更多，PACF则为AR检验服务的更多，这里非常不专业的归纳一下，若ACF截尾则倾向于MA现象显著，若PACF截尾则倾向于AR现象显著。
2.Ljung-Box.test检验
Ljung-Box倾向于检验总体平稳序列相关性，原假设：无自相关。
,
时拒绝原假设，其中
是自由度为k的卡方分布的
分位点
3.t.test
对平稳序列进行均值为0的t检验
这是为了mean equation中include.mean=TRUE/FALSE做判断。
1.2.2 初步定阶与参数估计[innovation is Gaussian]
1.2.2.1 初步定阶
初步定阶普遍采用eacf法，这个定阶是不唯一的，从而我们会得到多个model，但都要对他们进行参数估计，model有效性的检验与筛选在1.2.3节会具体说明。
1.2.2.2 参数估计
ARMA model参数估计这部分只粗略提一下方法，其中需要注意的唯一一点是，innovation是Gaussian分布，要用Gaussian PDF.
1.OLS/CSS 最小二乘
略
2.MLE 极大似然估计
- not conditional MLE
- exact MLE
- conditional MLE
3.CSS-ML
初值用OLS/CSS法确定，arima参数的计算用MLE
1.2.3模型筛选+模型诊断(model diagnostics)
该部分主要分为两大部分：模型比较与选择，模型残差诊断。
模型比较与选择(P92)
样本内比较
- AIC/BIC 信息准则 小者优
- 比较模型残差的方差估计 小者优
- 实际意义的比较
样本外比较
- 回测检验
将样本分为训练集和回测集，利用训练集所估计完参数的模型来对回测集进行预测，计算与真实值之间的均方误差(MSFE)
，或者平均绝对预测误差(MAFE)也OK。
模型残差诊断
即检验残差是否服从白噪声预设分布假设。
1.检验残差分布 - 残差序列图 - Q-Q plots - 利用相关分布检验方法，检验残差是否符合预设的Gaussian分布
2.检验残差自相关 - ACF - Ljung-Box test
1.2.4 得到最优ARMA model
1.3 GARCH效应的检验与建模
1.3.1 GARCH检验
对序列
进行Ljung-Box test/Lagrange Multiplier test的序列相关检验,原假设都为序列不相关，若拒绝则具有GARCH效应。
1.3.2 GARCH效应的参数估计与阶数选择
通常来说，GARCH效应一般不会太高阶，金融数据的GARCH效应大致为GARCH(1,0),(0,1),(1,1),(1,2),(2,1),(2,2)等。因此可以将这几个阶数分别代入GARCH模型中进行参数估计。
1.3.2.1 参数估计
一般采用极大似然估计进行参数估计（什么分布用什么PDF）
1.3.2.2 阶数选择
- AIC/BIC 信息准则
1.3.3模型诊断(model diagnostics)
为了诊断GARCH模型的有效性，对标准化后的残差
进行ACF,PACF,Ljung-Box检验自相关，然后再检验其分布是否符合预设分布。
1.4 求出k期均值与方差预测
我有一个疑问 对于处于在已知期(0,t)之间的
，存在下列等式嘛?还是说已知期的
是可以根据总体GARCH模型得到的？：
对于总体模型有：
1.4.1 均值预测
从而我们得到k期均值预测：
1.4.2 方差预测
由
的ARMA形式：
可得
的
形式:
从而可得第k期的预测误差forecast error
满足：
从而得到k期内的预测误差
满足：
由上面计算的预测误差，我们就可以计算出预测期的expected volatility:
所以接下来只需求出每期volatility的预测即可：
这里的
是通过对真实值
在
的条件下等式两边去期望所得。 后面的
亦是如此。 从而，我们求得了各期volatility的预测，从而也得到了k期内的volatility的预测。
1.5 Calculate the VaR
2 基于innovation服从标准学生t分布下的算法——仿真法(simulation)
2.1 前三步
前三步大致相同，但注意在ARMA和GARCH的参数估计和模型诊断时要记住innovation是服从标准学生t分布的。
2.2 仿真法(simulation)
2.2.1 第1期均值与方差预测
1.
2.随机产生一个标准t分布的扰动项
作为
的仿真值，从而得到仿真的新息innovation
3.再由新得到的仿真的innovation
,即可得到t+1期的loss variable仿真值：
2.2.2 第2期均值与方差预测
1.
2.再随即产生一个标准t分布的扰动项
作为对
的仿真值，从而得到t+2期的仿真的新息innovation
3.接着计算loss variable的仿真值：
2.2.3 接着计算第k期均值与方差预测,并得到k期内的均值与方差预测
2.2.4 计算第一次仿真的VaR
2.2.5 重复 2.2.1-2.2.4 仿真n次，得到n个VaR
重复n次仿真后，我们得到了
从而得到最终仿真结果：
3 CODE IN R
这里直接给出rmd文档代码。有相关注释，如果需要我可以私发PDF版
---
title: "Financial Econometrics Assignment4"
author: "41621069 MaoHua Wang"
date: "2018.12.23"
output:
html_document:
df_print: paged
editor_options:
chunk_output_type: console
---
# Exercise 1
Consider a long position of $1 million on the Apple stock. To assess the risk of the position, we employ daily returns of the stock from January 2, 2001 to September 30, 2011 for 2704 observations. The daily simple returns are obtained from CRSP and in the file d-aapl-0111.txt. Let the tail probability be p=0.01. Compute the VaR of the position for the next trading day and the next 10 trading days using the following methods:
(a) The RiskMetrics method. Write down the fitted special IGARCH(1,1) model.
(b) A Gaussian GARCH model. Write down the fitted model.
(c) A GARCH model with a standardized Student-t innovation. Write down the fitted model. [You should use simulation to compute VaR for the next 10 trading days based on the fitted model.]
```{r include=FALSE}
library(fGarch)
library(TSA)
library(rugarch)
library(fUnitRoots)
setwd('C:/Users/82335/Desktop/grade3/Financial Econometrics/assignment4')
data_1 = read.table('d-aapl-0111.txt',header = TRUE)
```
## (a) Risk Metrics
In Risk Metrics, the loss variable $a_t$ is logreturn with the Gaussian distribution(0,$sigma_t^2$) , and we assume that the volatility $sigma_tsim IGARCH(1,1)$,which means that:
$$mean quad equation: x_t = 0+a_t, a_t = sigma_tvarepsilon_t $$
$$volatility quad equation: sigma_t^2 = alphasigma_{t-1}^2 + (1-alpha)a_{t-1}^2, 0<alpha<1$$
```{r warning=FALSE}
#get the log return
data_1$log_return = log(data_1$rtn+1)
#Estimation in IGARCH(1,1)
mspec_a = ugarchspec(mean.model = list(armaOrder=c(0,0),include.mean=FALSE),variance.model = list(model='iGARCH',garchOrder=c(1,1)),fixed.pars = list(omega = 0))
fitmodel_a = ugarchfit(data = -data_1$log_return,spec = mspec_a)
#we get the alpha = 0.0345
#forecast
forcast_a1 = ugarchforecast(fitmodel_a,n.ahead = 1)
sigma_a1 = forcast_a1@forecast$sigmaFor[1]
VaR_a1 = 0 + 1*qnorm(0.99)*sigma_a1 #1 is mean 1 million
#Square root of time rule
VaR_a10 = sqrt(10)*VaR_a1
```
```{r}
VaR_a1
```
```{r}
VaR_a10
```
## (b)
### 1. Weak stationary test
```{r echo=TRUE, warning=FALSE}
ws = adfTest(data_1$log_return) #reject H0, proof the log_return series are weak stationary.
ws@test$p.value
```
### 2. ARMA model
```{r echo=FALSE, fig.height=3, fig.width=9, message=FALSE, warning=FALSE}
acf(data_1$log_return,main='acf')
pacf(data_1$log_return,main='pacf')
```
```{r}
box = Box.test(data_1$log_return,lag=12,type='Ljung-Box')
box$p.value
```
reject H0, exist series correlation
```{r}
t = t.test(data_1$log_return)
t$p.value
```
reject H0, exist mean, so the model need 'include.mean=TRUE'
```{r}
eacf(data_1$log_return) #ARMA(0,4)
```
### 3. GARCH model
```{r warning=FALSE}
#calculate the resid^2
for (i in 1:length(data_1$log_return)){
data_1$resid2[i] = (data_1$log_return[i]-mean(data_1$log_return[1:i]))^2
}
box = Box.test(data_1$resid2)#reject H0, i.e. exist GARCH effect
#choose GARCH order
order_b = c(c(1,0),c(0,1),c(1,1),c(1,2),c(2,1),c(2,2))
for (i in 1:6){
mspec_b = ugarchspec(mean.model = list(armaOrder=c(0,4),include.mean=TRUE),variance.model = list(garchOrder=c(order_b[2*i-1],order_b[2*i])))
fit_b = ugarchfit(data=data_1$log_return,spec=mspec_b)
aic_b = -2*fit_b@fit[["LLH"]]+2*length(fit_b@fit[["solver"]][["sol"]][["pars"]])
}
#GARCH(1,1) is the best
```
### 4. Forecast
Volatility forecast of the h-period loss variable:
$$begin{split}
V_t(e_t[k])&=V_t[sum_{i=1}^ke_t(i)]
&=sum_{i=1}^kV_t[e_t(i)]
&=quad(1+psi_{1}+psi_{2}+cdots+psi_{k-1})^2V_t(a_{t+1})
&quad+(1+psi_{1}+psi_{2}+cdots+psi_{k-2})^2V_t(a_{t+2})
&quad+cdots
&quad+(1+psi_1)^2V_t(a_{t+k-1})
&quad+V_t(a_{t+k})
&=quad(1+psi_{1}+psi_{2}+cdots+psi_{k-1})^2sigma^2_t(1)
&quad+(1+psi_{1}+psi_{2}+cdots+psi_{k-2})^2sigma^2_t(2)
&quad+cdots
&quad+(1+psi_1)^2sigma_t^2(k-1)
&quad+sigma_t^2(k)
end{split}$$
Where $psi_i$ is the parameter in $x_tsim MA(infty)$.
i.e.
$$V[h]=(sum_{i=0}^{h-1} psi_i)^2 sigma^2_1+....+(1+psi_1)^2 sigma^2_{h-1}+sigma^2_h$$
```{r include=FALSE}
i = 3
mspec_b = ugarchspec(mean.model = list(armaOrder=c(0,4),include.mean=TRUE),variance.model = list(garchOrder=c(order_b[2*i-1],order_b[2*i])))
fit_b = ugarchfit(data=data_1$log_return,spec=mspec_b)
forecast_b=ugarchforecast(fit_b,n.ahead=10)
psi = c(1:10)
psi[1] = 1
psi[6:10] = 0
for (i in 2:5){
psi[i] = fit_b@fit$coef[i]
}
VaR_b1 = 1*(forecast_b@forecast$seriesFor[1]+qnorm(0.99)*forecast_b@forecast$sigmaFor[1])
sigma2_10 = 0
for (i in 1:10){
sum_k = 0
for (j in 1:i){
sum_k = sum_k+psi[j]
}
sigma2_10 = sigma2_10+(sum_k)^2*(forecast_b@forecast$sigmaFor[11-i])^2
}
loss_variable_b = sum(forecast_b@forecast$seriesFor)
sigma_10 = sqrt(sigma2_10)
VaR_b10 = 1*(loss_variable_b+qnorm(0.99)*sigma_10)
```
```{r}
VaR_b1
VaR_b10
```
## (c)
As we know, the first 3 steps is same as above in the GARCH model with a standardized Student-t innovation. So we just need change the distribution of the white noise when setting.
Here's the steps:
### Simulate 1-step-ahead
1.$GARCH to sigma^{ast2}_{t+1}=alpha_0 + sum_{i=1}^ualpha_ia_{t+1-i}^2+sum_{j=1}^vbeta_jsigma_{t+1-j}^2tosigma^{ast2}_{t+1}$
2.Simulate a noise which obeys normal distribution: $varepsilon^{ast}_{t+1}$ as$varepsilon_{t+1}$ then we get the simulated innovation $a^{ast}_{t+1}=sigma^{ast}_{t+1}varepsilon^{ast}_{t+1}$
3.By the simulated innovation$a^{ast}_{t+1}$, we can get the simulation of loss variable at t+1: $$x_{t+1}^{ast} = phi_o+sum_{i=1}^pphi_ix_{t+1-i}+a^{ast}_{t+1}+sum_{j=1}^qtheta_ja_{t+1-j}$$
### Simulate 2-step-ahead
1.$sigma_{t+2}^{ast2}=alpha_0+alpha_1a_{t+1}^{ast2}+sum_{i=2}^ualpha_ia_{t-i+2}^{2}+beta_1sigma_{t+1}^{ast2}+sum_{j=2}^vbeta_jsigma_{t-j+2}^2$
2.$a^{ast}_{t+2}=sigma^{ast}_{t+2}varepsilon^{ast}_{t+2}$
3.$x_{t+2}^{ast} = phi_o+phi_1x_{t+1}^{ast}+sum_{i=2}^pphi_ix_{t+2-i}+a^{ast}_{t+2}+theta_1 a^{ast}_{t+1}+sum_{j=2}^qtheta_ja_{t+2-j}$
### Simulate k-step-ahead
$x_t[k] = x_{t+1}^{ast}+x_{t+2}^{ast}+cdots+x_{t+k}^{ast}$
$V_t(e_t[k])=sigma_{t+k}^{ast2}+(1+psi_1)^2sigma_{t+k-1}^{ast2}+cdots+(1+psi_1+cdots+psi_{k-1})^2sigma_{t+1}^{ast2}$
### Calculate the VaR at the first time
$VaR^{(1)}_{1-p} = x_t[k]+t_{1-p}^{ast}V_t(e_t[k])$
### Repeat n times(such as 10000)
cycle-cycle-cycle...
## Code in R
```{r}
#build model
mspec_c=ugarchspec(mean.model=list(armaOrder=c(0,4),include.mean=TRUE),variance.model=list(garchOrder=c(1,1)),distribution.model = "std")
fit_c=ugarchfit(data=data_1$log_return,spec=mspec_c)
psi = c(1:10)
psi[1] = 1
psi[6:10] = 0
for (i in 2:5){
psi[i] = fit_c@fit$coef[i]
}
#get t df
df=fit_c@fit$coef[9]
#get a_t and ahead a_{t-1},a_{t-2},a_{t-3}
resid = as.numeric(residuals(fit_c))
VaR = matrix(0,10000,1)
for (t in 1:10000){
a_t_0 = tail(resid,4)[4]
a_t_1 = tail(resid,4)[3]
a_t_2 = tail(resid,4)[2]
a_t_3 = tail(resid,4)[1]
a_t = c(1:14)
a_t[1] = a_t_3
a_t[2] = a_t_2
a_t[3] = a_t_1
a_t[4] = a_t_0
h=10
#set forecast variance set
forecast_variance=matrix(0,h,1)
forecast_xt=matrix(0,h,1)
sigmahat=as.numeric(sigma(fit_c))
sigma.last.period=tail(as.numeric(sigma(fit_c)),1)
fit_c@fit$coef
for (i in 1:h){
newsigma2 = fit_c@fit$coef[6]+fit_c@fit$coef[7]*a_t[i+3]^2+fit_c@fit$coef[8]*sigma.last.period^2
forecast_variance[i] = newsigma2
a_t[i+4] = sqrt(as.numeric(newsigma2))*rstd(1,mean=0,sd=1,nu=df)
sigma.last.period = sqrt(newsigma2)
forecast_xt[i] = fit_c@fit$coef[1]+a_t[i+4]+fit_c@fit$coef[2]*a_t[i+3]+fit_c@fit$coef[3]*a_t[i+2]+fit_c@fit$coef[4]*a_t[i+1]+fit_c@fit$coef[5]*a_t[i]
}
variance10 = 0
for (i in 1:10){
sum_k = 0
for (j in 1:i){
sum_k = sum_k+psi[j]
}
variance10 = variance10+(sum_k)^2*forecast_variance[11-i]
}
xt10=sum(forecast_xt)
sigma10 = sqrt(variance10)
VaR[t] = 1*(xt10+qt(0.99,df)*sigma10)
}
mean(VaR)
```
# Exercise 2
Cointegration analysis: the file q-4macro.txt contains two U.S. quarterly macroeconomic variables and two interest rates. Construct a multivariate time series $z_t$, where $z_{1t}$ is the logarithm of the U.S. real gross national product (lngnp), $z_{2t}$ is the rate of interest rate of U.S. 3-month treasury bills (tb3m), $z_{3t}$ is logarithm of the U.S. M1 money stock (lnm1), and $z_{4t}$ is interest rate of U.S. 10-year government security (gs10).
(a) Are these four series unit-root stationary? Why?
(b) What is the optimal VAR model for $z_t$ by sequential likelihood ratio tests, AIC, BIC, and HQ, respectively?
(c) Please identify any cointegration among $z_t$. If there cointegration exists, how many cointegrating vectors of $z_t$ and what are they?
(d) Please estimate the Error Correction Model (VECM) of $z_t$ and write down the fitted VEM model.
## (a)
```{r echo=FALSE, fig.height=3, fig.width=9, message=FALSE, warning=FALSE, paged.print=FALSE}
data_2 = read.table('q-4macro.txt',header = TRUE)
zt=cbind(data_2$lngnp,data_2$tb3m,data_2$lnm1,data_2$gs10)
colnames(zt)=c("lngnp","tb3m","lnm1","gs10")
plot(zt[,1],type='l',main='lngnp')
#adfTest(zt[,1],type = 'c')@test$p.value
plot(zt[,2],type='l',main='tb3m')
#adfTest(zt[,1],type = 'nc')@test$p.value
plot(zt[,3],type='l',main='lnm1')
#adfTest(zt[,1],type = 'c')@test$p.value
plot(zt[,4],type='l',main='gs10')
#adfTest(zt[,1],type = 'nc')@test$p.value
```
Obviously, these four series are not unit-root stationary.
## (b)
```{r}
library(MTS)
m_b = VARorder(zt,12)
```
AS we can see from the VARorder plot:
$$begin{array}{c|ccc}
Information Criterion & text{Order Select}
hline
AIC & 6
BIC & 2
HQ & 2
p-value & 11
end{array}$$
Suppose we use the AIC select 6.
## (c)
```{r include=FALSE}
# test cointegration
library(urca)
```
```{r}
# trace test(type='trace')
# pairs trade with consistent(ecdet='const')
# VAR has 2 subject
fit_c=ca.jo(zt,type="trace",ecdet="const",K=2,spec="transitory")
summary(fit_c)
```
The test.value of r = 0 is 122.19>60.16, reject H0.
The test.value of r = 1 is 40.82>34.91, reject H0.
The test.value of r = 2 is 13.92<17.85, accept H0.
So there exists 2 cointegrated vector:
$$beta = begin{bmatrix} 1.0000 &1.0000 -0.4091&-0.6139 -0.8585 &-0.6857 0.4316&0.5926 end{bmatrix}$$
## (d)
```{r}
vecm_d=cajorls(fit_c,r=2) # assume rank(pi)=m=2
show(vecm_d)
```
From the ECM equation:
$$Delta r_t = mu_t +alphabeta'r_{t-1}+phi_1^{ast}Delta r_{t-1}+cdots+phi_{p-1}^{ast}Delta r_{t-p-1}+varepsilon_t$$
We can get the model below:
$$
begin{bmatrix}
Delta lngnp
Delta tb3m
Delta lnm1
Delta gs10
end{bmatrix}
=
begin{bmatrix}
0.0056 &-0.0022
0.0644&-0.0601
-0.0042 &-0.0015
-0.2384&0.1134
end{bmatrix}
beta'
begin{bmatrix}
lngnp_{t-1}
tb3m_{t-1}
lnm1_{t-1}
gs10_{t-1}
end{bmatrix}
+
begin{bmatrix}
0.2047&0.0033&-0.0430&-0.0025
11.68887&0.1820&5.2033&0.0876
-0.0492&-0.0021&0.5759&-0.0043
6.7372&-0.0554&6.9089&0.2381
end{bmatrix}
begin{bmatrix}
Delta lngnp_{t-1}
Delta tb3m_{t-1}
Delta lnm1_{t-1}
Delta gs10_{t-1}
end{bmatrix}
+
begin{bmatrix}
varepsilon_{1t}
varepsilon_{2t}
varepsilon_{3t}
varepsilon_{4t}
end{bmatrix}
$$
$$where quad beta'=begin{bmatrix}
1.0000 &-0.4091&-0.8585&0.4316
1.0000&-0.6139&-0.6857&0.5926
end{bmatrix}$$
3.1 TEST
3.1.1 ADF test
library(fUnitRoots)
adfTest(ln.rtn,lags=12,type=c("nc"))
type=c(): choose the series type: random walk, random walk with drift, trend-stationary time series.The default is 'c'.
1. nc: random walk. no intercept(constant) nor time trend.
2. c: random walk with drift.with an intercept(constant) but no time trend
3. ct: trend stationary time series. with an intercept(constant) and a time trend.
3.1.2 Ljung-Box test
3.1.3 Lagrange Multiplier test
3.1.4 t test
3.1.5 Gaussian distribution test
3.1.6 standard t distribution test
3.2 function
3.2.1 ACF
3.2.2 PACF
3.2.3 EACF
3.3 MODEL
3.3.1 ARMA
3.3.2 GARCH
3.4 Algorithm
3.4.1 Risk Metrics
3.4.2 Econometric Approach
3.4.3 Quantile Estimation
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
