Analysis of Financial Time Series in Stan - Chapter 10 - Multivariate Volatility Models and Their Applications
In this series we are trying to reproduce the models and examples listed in the book “Analysis of Financial Time Series, 3rd Edition”, by Ruey S. Tsay, using Stan https://mc-stan.org/ and the package RStan https://cran.r-project.org/web/packages/rstan/index.html. The main repository for the presented models and data can be found at https://github.com/marcomarconi/AFTS_with_Stan.
{
require(tidyverse)
require(rstan)
rstan_options(auto_write = TRUE)
rstan_options(javascript = FALSE)
theme_set(theme_classic(base_size = 24))
}10.1 EXPONENTIALLY WEIGHTED ESTIMATE
The exponentially weighted moving-average (EWMA) is way to estimate the unconditional covariance matrix of the innovations of a time series. For a sufficiently long time series it can be approximated as
\[ \widehat{\Sigma}_{t} = (1-\lambda)a_{t-1}a'_{t-1}+\lambda\widehat{\Sigma}_{t-1} \] This model can be written in Stan as:
data {
int<lower=0> N;
int<lower=0> K;
vector[K] y[N];
matrix[K,K] Sigma1;
int approx; // use the approximated formula (if N is large)
}
parameters {
real<lower=0,upper=1> lambda;
vector[K] mu;
}
transformed parameters {
matrix[K,K] Sigma[N];
Sigma[1] = Sigma1;
if(approx) { // approximated formula (if N is large)
for (t in 2:N)
Sigma[t] = (1-lambda) * ((y[t-1]-mu) * to_row_vector(y[t-1]-mu)) + lambda * Sigma[t-1];
}else{ // recursive formula
for (t in 2:N) {
Sigma[t] = Sigma1 ;
for(j in 1:(t-1))
Sigma[t] = Sigma[t] + pow(lambda,(j-1)) * ((y[t-1]-mu) * to_row_vector((y[t-1]-mu))) ;
Sigma[t] = ((1-lambda) / (1-pow(lambda,(t-1)))) * Sigma[t];
}
}
}
model {
lambda ~ beta(1, 1);
for (t in 2:N)
y[t] ~ multi_normal(mu, Sigma[t]);
}
Example 10.1.
We are going to fit a univariate GARCH model to daily log returns of the Hang Seng index of Hong Kong and the Nikkei 225 index of Japan from January 4, 2006, to December 30, 2008, for 713 observations, and later smooth the estimated volatility using EWMA.
x <- read.table("data/d-hkjp0608.txt", header=T)
hk <- x$HK %>% log %>% diff %>% `*`(100)
jp <- x$JP %>% log %>% diff %>% `*`(100)
r <- cbind(hk, jp)
plot.ts(r, main="log returns", cex.axis=2, cex.lab=2, cex.main=2)
Let’s see the GARCH estimate first:
GARCH <- stan_model("models/Chapter_3/GARCH.stan")
fit_HK <- sampling(GARCH, data = list(N = length(hk), y=hk, sigma1=sd(hk), K=1), chains = 1, cores = 4, refresh=0)
fit_JP <- sampling(GARCH, data = list(N = length(jp), y=jp, sigma1=sd(jp), K=1), chains = 1, cores = 4, refresh=0)
pars1 <- rstan::extract(fit_HK, pars="sigma")[[1]]
pars2 <- rstan::extract(fit_JP, pars="sigma")[[1]]
plot.ts(cbind(hk=colMeans(pars1), jp=colMeans(pars2)), main="", cex.axis=2, cex.lab=2, )
and now apply the EWMA smoothing to the estimated volatilities, you can see that the parameters are similar to Tsay’s:
EWMA <- stan_model("models/Chapter_10/EWMA.stan")
fit_EWMA <- sampling(EWMA, data = list(N = nrow(r), y=r, Sigma1=cov(r), K=2, approx=1), chains = 4, cores = 4)EWMA <- stan_model("models/Chapter_10/EWMA.stan")
fit_EWMA <- readRDS("models/Chapter_10/saved_RDS/fit_EWMA.rds")print(fit_EWMA, pars=names(fit_EWMA)[1:3])## Inference for Stan model: EWMA.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## lambda 0.93 0 0.01 0.92 0.93 0.93 0.93 0.94 2864 1
## mu[1] 0.07 0 0.05 -0.02 0.04 0.07 0.11 0.16 2600 1
## mu[2] -0.01 0 0.05 -0.11 -0.04 -0.01 0.02 0.08 2437 1
##
## Samples were drawn using NUTS(diag_e) at Sun Jul 18 12:50:08 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Sigma <- rstan::extract(fit_EWMA, pars="Sigma")[[1]]
Sigma <- colMeans(Sigma) %>% apply(c(1,2), mean)
colnames(Sigma) <- c("hk", "jp")
plot.ts(Sigma, main="", cex.axis=2, cex.lab=2, )
10.2.1 Diagonal Vectorization (VEC) Model
The DVEC model is a generalization of the EWMA proposed by Bollerslev, Engle, and Wooldridge (1988). It can be coded in Stan as follows:
data {
int<lower=0> N;
int K;
vector[K] y[N];
matrix[K,K] Sigma1;
}
parameters {
vector[K] mu;
real<lower=0> A11;
real<lower=0> A21;
real<lower=0> A22;
real<lower=0> arch11;
real<lower=0> arch21;
real<lower=0> arch22;
real<lower=0> garch11;
real<lower=0> garch21;
real<lower=0> garch22;
}
transformed parameters {
matrix[K,K] Sigma[N];
Sigma[1] = Sigma1;
for(t in 2:N) {
vector[K] a = y[t-1] - mu;
Sigma[t][1,1] = A11 + arch11 * a[1]*a[1] + garch11 * Sigma[t-1][1,1];
Sigma[t][2,1] = A21 + arch21 * a[1]*a[2] + garch21 * Sigma[t-1][2,1];
Sigma[t][1,2] = Sigma[t][2,1];
Sigma[t][2,2] = A22 + arch22 * a[2]*a[2] + garch22 * Sigma[t-1][2,2];
}
}
model {
mu ~ normal(0,1);
[A11, A21, A22] ~ normal(0,0.001);
[arch11, arch21, arch22, garch11, garch21, garch22] ~ normal(0,1);
for(t in 2:N)
y[t] ~ multi_normal(mu, Sigma[t]);
}
Fitting this model is quite hard as the positive-definitiveness of covariance matrix may not be guaranteed.
Let’s fit it to the monthly returns of Pfizer and Merck stocks:
r <- read.table("data/m-pfemrk6508.txt", header=T)
plot.ts(r[,-1], main="", cex.axis=2, cex.lab=2)
DVEC <- stan_model("models/Chapter_10/DVEC.stan")
fit_DVEC <- sampling(DVEC, data = list(N = nrow(r), y=r[,-1], Sigma1=cov(r[,-1]), K=2), iter=500,chains = 1, cores = 4, control=list(adapt_delta=0.99))We can plot the estimated volatilities and time-varying correlations. Unfortunately, the time-varying correlation do not match the one presented in the book:
Sigma <- rstan::extract(fit_DVEC, pars="Sigma")[[1]]
sigmas <- colMeans(Sigma) %>% apply(c(1,2), mean)
colnames(sigmas) <- c("pfe", "mrk")
plot.ts(sigmas, main="", cex.axis=2, cex.lab=2, )
corr <- apply(Sigma, c(2,3,4), mean) %>% apply(.,1,cov2cor) %>% t
plot.ts(corr[,2])
10.4 GARCH MODELS FOR BIVARIATE RETURNS
Univariate volatility models can be readily extended to the multivariate case. In this situation, the conditional covariatence matrix at a given time can be reparametrized as :
\[ \Sigma_t = D_t\rho_tD_t \]
where \(\rho_t\) is the conditional correlation matrix of \(a_t\) , and \(D_t\) is a k × k diagonal matrix consisting of the conditional standard deviations of elements of \(a_t\). Therefore, to model the volatility is suffices to consider the k(k + 1)/2-dimensional vector:
\[ \Xi_{t} = (\sigma_{11,t},...,\sigma_{kk,t},\zeta'_t) \]
where \(\zeta'_t\) is a k(k − 1)/2-dimensional vector obtained by stacking columns of the correlation matrix \(\rho_t\) , but using only elements below the main diagonal.
10.4.1 Constant-Correlation Models
Under the assumption of constant correlation the model the volatility model becomes:
\[ \Xi_t = (\sigma_{11,t},\sigma_{22,t})' \\ \Xi_t = \alpha_0+\alpha_1a^2_{t-1}+beta_1\Xi_{t-1} \]
A simple bivariate M-GARCH model of this kind can be modelled in Stan as:
data {
int<lower=0> N;
vector[2] y[N];
matrix[2,2] Sigma1;
}
parameters {
vector[2] mu;
real<lower=-1,upper=1> rho;
vector<lower=0>[2] alpha0;
cov_matrix[2] Alpha1;
cov_matrix[2] Beta1;
}
transformed parameters {
vector[2] Xi[N];
Xi[1][1] = Sigma1[1,1];
Xi[1][2] = Sigma1[2,2];
for (t in 2:N) {
vector[2] a_1 = y[t-1] - mu;
Xi[t] = alpha0
+ Alpha1 * [a_1[1]*a_1[1], a_1[2]*a_1[2]]'
+ Beta1 * Xi[t-1];
}
}
model {
matrix[2,2] Rho ;
mu ~ normal(0, 1);
alpha0 ~ normal(0, 1);
Rho = [[1, rho], [rho, 1]];
for (t in 2:N)
y[t] ~ multi_normal(mu, quad_form_diag(Rho, Xi[t]));
}
generated quantities {
vector[N] log_lik;
matrix[2,2] Rho = [[1, rho], [rho, 1]];
for (t in 1:N)
log_lik[t] = multi_normal_lpdf(y[t] | mu, quad_form_diag(Rho, Xi[t]));
}Let’s fit this model to the daily log returns of Hong Kong and Japanese markets seen before:
x <- read.table("data/d-hkjp0608.txt", header=T)
hk <- x$HK %>% log %>% diff %>% `*`(100)
jp <- x$JP %>% log %>% diff %>% `*`(100)
r <- cbind(hk, jp)fit_MGARCH_hkjp <- stan_model("models/Chapter_10/fit_MGARCH_hkjp.stan")
fit_MGARCH_hkjp <- sampling(MGARCH, data = list(N= nrow(r), y=r, Sigma1=cov(r)), iter=500,chains = 4, cores = 4, init=0)The parameters estimated by the model are close to Tsay’s:
print(fit_MGARCH_hkjp, pars=names(fit_MGARCH_hkjp)[1:13])## Inference for Stan model: MGARCH.
## 1 chains, each with iter=500; warmup=250; thin=1;
## post-warmup draws per chain=250, total post-warmup draws=250.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 0.10 0 0.06 0.00 0.06 0.10 0.13 0.21 128 1
## mu[2] 0.00 0 0.05 -0.08 -0.03 0.00 0.04 0.10 141 1
## rho 0.67 0 0.02 0.63 0.65 0.67 0.68 0.70 326 1
## alpha0[1] 0.22 0 0.03 0.17 0.20 0.22 0.24 0.30 231 1
## alpha0[2] 0.14 0 0.03 0.09 0.11 0.13 0.15 0.20 157 1
## Alpha1[1,1] 0.05 0 0.01 0.04 0.04 0.05 0.05 0.06 204 1
## Alpha1[2,1] 0.00 0 0.00 -0.01 0.00 0.00 0.01 0.01 153 1
## Alpha1[1,2] 0.00 0 0.00 -0.01 0.00 0.00 0.01 0.01 153 1
## Alpha1[2,2] 0.02 0 0.01 0.01 0.02 0.02 0.03 0.04 153 1
## Beta1[1,1] 0.75 0 0.03 0.69 0.73 0.75 0.78 0.81 153 1
## Beta1[2,1] 0.02 0 0.03 -0.03 0.00 0.02 0.04 0.07 115 1
## Beta1[1,2] 0.02 0 0.03 -0.03 0.00 0.02 0.04 0.07 115 1
## Beta1[2,2] 0.84 0 0.03 0.77 0.81 0.84 0.86 0.90 105 1
##
## Samples were drawn using NUTS(diag_e) at Sat Jul 24 16:30:02 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).pars <- rstan::extract(fit_MGARCH_hkjp)
sigmas <- pars$Xi %>% colMeans() %>% sqrt
colnames(sigmas) <- c("hk", "jp")
plot.ts(sigmas, cex.axis=2, cex.lab=2)
The Ljung–Box statistics of the standardized residuals do not show any sort of model inadequancy:
library(portes)
mu <- pars$mu %>% colMeans()
a <- r - mu
residuals <- a / sigmas
LjungBox(residuals, lags = 1:12)## lags statistic df p-value
## 1 16.67645 4 0.002233669
## 2 19.13570 8 0.014150553
## 3 24.15928 12 0.019349593
## 4 25.11302 16 0.067857624
## 5 28.67279 20 0.094404890
## 6 29.36183 24 0.206810302
## 7 34.40102 28 0.187980357
## 8 40.01167 32 0.156211919
## 9 44.12083 36 0.165878608
## 10 53.05252 40 0.081044748
## 11 54.30802 44 0.137200146
## 12 58.63381 48 0.13982874210.4.2 Time-Varying Correlation Models
Constant correlation is quite a strong assumption, especially for long term relationships like between IBM stock and the S&P 500 index. We are going to apply a time-varying correlation model that use the Cholesky decomposition of the conditional covariance matrix. The formulas are quite involved so please refer to equations (10.28) and (10.30) of the book. The model can be coded in Stan as follows, notice that the model accept three type of likelihoods: 1 - using a multivariate normal with the covariance decomposed as LGL’ Cholesky, 2 - using the special multivariate normal for Cholesky provided by Stan, and 3 - using the analytic likelihood described in the book at equqation (10.20). We are going to use the second option in this example:
data {
int<lower=0> N;
vector[2] y[N];
vector[2+1] Xi0;
int likelihood;
int prior;
}
parameters {
vector[2] mu;
real<lower=0> alpha10;
real<lower=0> alpha11;
real<lower=0> alpha21;
real<lower=0> alpha22;
real<lower=0> alpha20;
real<lower=0> beta11;
real<lower=0> beta21;
real<lower=0> beta22;
real lambda0;
real lambda1;
real lambda2;
}
transformed parameters {
vector[2+1] Xi[N];
vector[2] a[N];
vector[2] b[N];
Xi[1,] = Xi0;
a[1] = y[1] - mu;
b[1][1] = a[1][1];
b[1][2] = a[1][2] - Xi[1,3]*a[1][1] ;
// equation (10.28)
for (t in 2:N) {
Xi[t,1] = (alpha10 + alpha11 * pow(b[t-1][1], 2) + beta11 * Xi[t-1,1]); // g11
Xi[t,2] = (alpha20 + alpha21 * pow(b[t-1][1], 2) + alpha22 * pow(b[t-1][2], 2) + beta21 * Xi[t-1,1] + beta22 * Xi[t-1,2]); // g22
Xi[t,3] = lambda0 + lambda1 * Xi[t-1,3] + lambda2 * a[t-1][2]; // q21
a[t] = y[t] - mu;
b[t][1] = a[t][1];
b[t][2] = a[t][2] - Xi[t,3]*a[t][1] ;
}
}
model {
matrix[2,2] G;
matrix[2,2] L = diag_matrix(rep_vector(1,2));
matrix [2,2] Sigma;
mu ~ normal(0, 2);
[lambda0, lambda1, lambda2] ~ normal(0, 1);
[alpha10, alpha11, alpha20, alpha21, alpha22] ~ normal(0, 2);
[beta11, beta21, beta22] ~ normal(0, 1);
for (t in 2:N) {
G[1,1] = (Xi[t,1]);
G[2,2] = (Xi[t,2]);
G[2,1] = 0;
G[1,2] = 0;
L[2,1] = Xi[t,3];
if(likelihood == 1) // multi_normal
Sigma = L * G * L';
else if(likelihood == 2) // multi_normal_cholesky
Sigma = L * sqrt(G);
if(!prior){
if(likelihood == 1) // multi_normal
y[t] ~ multi_normal(mu, Sigma);
else if(likelihood == 2) // multi_normal_cholesky
y[t] ~ multi_normal_cholesky(mu, Sigma);
else if(likelihood == 0) { // book's loglik
real tmp = 0;
for(k in 1:2)
tmp += log(Xi[t][k]) + pow(b[t][k],2) / Xi[t][k];
target += -0.5 * tmp;
}
}
}
}
generated quantities {
vector[N] rho;
vector[N] log_lik;
matrix[2,2] G;
matrix[2,2] L = diag_matrix(rep_vector(1,2));
matrix [2,2] Sigma;
rho[1]=0.5;
for (t in 1:N) {
rho[t] = (Xi[t,3] * sqrt((Xi[t,1]))) / sqrt((Xi[t,2]) + pow(Xi[t,3],2) * (Xi[t,1])); // equation (10.30)
G[1,1] = (Xi[t,1]);
G[2,2] = (Xi[t,2]);
G[2,1] = 0;
G[1,2] = 0;
L[2,1] = Xi[t,3];
// we always use multi_normal_cholesky here
Sigma = L * sqrt(G);
log_lik[t] = multi_normal_cholesky_lpdf(y[t] | mu, Sigma);
}
}Let’s fit it to the S&P500 vs IBM data:
x <- read.table("data/m-ibmsp2608.txt", header=T)
r <- x[,2:3]*100fit_MGARCH_chol <- stan_model("models/Chapter_10/fit_MGARCH_chol.stan")
fit_MGARCH_chol <- sampling(MGARCH_chol, data = list(N= nrow(r), y=r, Xi0=c(sd(r[,1]), sd(r[,2]), cor(r[,1], r[,2])), prior=0, likelihood=2), iter=500,chains = 1, init=0)The parameters estimated are very close to Tsay’s
print(fit_MGARCH_chol, names(fit_MGARCH_chol)[1:13])## Inference for Stan model: MGARCH_cholesky.
## 1 chains, each with iter=500; warmup=250; thin=1;
## post-warmup draws per chain=250, total post-warmup draws=250.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu[1] 1.35 0.01 0.23 0.93 1.22 1.34 1.50 1.82 258 1.00
## mu[2] 0.61 0.01 0.15 0.35 0.51 0.62 0.71 0.89 246 1.00
## alpha10 3.75 0.08 0.89 2.09 3.14 3.76 4.26 5.57 130 1.01
## alpha11 0.12 0.00 0.03 0.07 0.10 0.12 0.14 0.18 117 1.01
## alpha21 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 90 1.00
## alpha22 0.09 0.00 0.02 0.06 0.08 0.09 0.10 0.12 141 1.00
## alpha20 0.29 0.01 0.13 0.10 0.19 0.27 0.35 0.61 238 1.00
## beta11 0.81 0.00 0.04 0.74 0.78 0.81 0.83 0.87 107 1.01
## beta21 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 275 1.00
## beta22 0.89 0.00 0.02 0.85 0.87 0.89 0.90 0.92 111 1.00
## lambda0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 177 1.01
## lambda1 0.99 0.00 0.00 0.99 0.99 0.99 0.99 1.00 142 1.01
## lambda2 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 183 1.01
##
## Samples were drawn using NUTS(diag_e) at Sun Jul 25 18:28:16 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Additionally, we can plot the fitted conditional correlation coefficient, with the 90% credible intervals:
library(lubridate)
rho <- rstan::extract(fit_MGARCH_chol, pars="phi")[[1]]
q <- t(apply(rho, 2, function(x) quantile(x, probs=c(0.05, 0.5, 0.95))))
q <- data.frame(Date=as_date(x[,1] %>% as.character()), q)
ggplot(q) + geom_ribbon(aes(Date, ymin=`X5.`, ymax=`X95.`), color="lightslateblue", fill="lightblue", size=0.2) + geom_line(aes(Date, `X50.`)) + ylab("Time-varying correlation") + xlab("") + ylim(c(0,1)) + theme(text=element_text(size=32))
As you can see, the estimated time-varying correlation is close to the one obtained in the book.
10.4.3 Dynamic Correlation Models
Tsay (2006) proposed a Dynamic Correlation Models where standardized innovations are assumed to follow a multivariate Student-t distribution and the diagonal matrix of the stardard deviations of the decomposed \(\Sigma\) seen before is defined as:
\[ D^2_t = \Lambda_0 + \Lambda_1D^2_{t-1}+\Lambda_2A^2_{t-1} \]
The equation can also be extended with a leverage effect element.
The correlation matrix is instead defined as:
\[
\rho_t = (1-\theta_1-\theta_2)\hat\rho+\theta_1\psi_{t-1}+\theta_2\rho_{t-1}
\] Where \(\hat\rho\) is the sample correlation matrix and \(\psi_t\) is the rolling sample correlation matrix of length \(m\).
This model is more parsimonious than the ones presented before as it uses the sample data in the specifications.
The DCC in Stan can be implemented as follows:
data {
int<lower=0> N;
int K;
vector[K] y[N];
matrix[K,K] Rho_hat;
matrix[K,K] Psi[N];
int prior;
}
parameters {
real<lower=0> nu;
vector<lower=1e-12>[K] alpha0;
vector<lower=0,upper=1>[K] alpha1;
vector<lower=0,upper=1>[K] alpha2;
real<lower=0,upper=1> theta1;
real<lower=0,upper=(1-theta1)> theta2;
}
transformed parameters {
matrix[K,K] D[N];
matrix[K,K] Rho_m[N];
matrix[N,K] a;
matrix[K,K] Alpha0 = diag_matrix(alpha0);
matrix[K,K] Alpha1 = diag_matrix(alpha1);
matrix[K,K] Alpha2 = diag_matrix(alpha2);
D[1] = diag_matrix([1, 1, 1, 1]');
Rho_m[1] = diag_matrix([1, 1, 1, 1]');
a[1,] = (y[1])';
for (t in 2:N) {
a[t,] = (y[t])';
D[t] = diag_matrix([1, 1, 1, 1]');
for(k in 1:K)
D[t,k,k] = sqrt(alpha0[k] + alpha1[k] * D[t-1,k,k]^2 + alpha2[k] * (a[t-1,k] * a[t-1,k]));
Rho_m[t] = (1-theta2-theta1) * Rho_hat + theta1 * Psi[t-1] + theta2 * Rho_m[t-1];
}
}
model {
matrix [K,K] Sigma;
nu ~ normal(0, 10);
alpha0 ~ normal(0,.1);
alpha1 ~ normal(0,1);
alpha2 ~ normal(0,.1);
[theta1, theta2] ~ normal(0, 1);
for (t in 1:N) {
Sigma = D[t] * Rho_m[t] * D[t];
if(!prior)
y[t] ~ multi_student_t(nu, rep_vector(0, K), Sigma);
}
}
generated quantities {
vector[N] log_lik;
for (t in 1:N) {
matrix [K,K] Sigma = D[t] * Rho_m[t] * D[t];
log_lik[t] = multi_student_t_lpdf(y[t] | nu, rep_vector(0, K), Sigma);
}
}We are going to fit the model to the daily exchange rates between U.S. dollar versus European euro and Japanese yen and the stock prices of IBM and Dell from January 1999 to December 2004:
x <- read.table("data/m-ibmsp2608.txt", header=T)
r <- x[,2:3]*100The model requires the rolling matrix correlations as input:
x <- read.table("data/d-fxsk9904.txt")*100
x <- as.matrix(x); colnames(x) <- c("USEU", "USJP", "IBM", "Dell")
N<-nrow(x)
K <- 4
m <- 69
Psi <- array(dim = c(N,K,K))
Psi[1,,] <- diag(rep(1, K))
a <- matrix(nrow=N, ncol=K)
a[1,] <- rep(0, K)
for(t in 2:N) {
a[t,] <- x[t,]
if(t > m)
Psi[t,,] <- cor(a[(t-m):(t-1),])
else
Psi[t,,] <- diag(rep(1, K))
}Let’s fit the model:
DCC <- stan_model("models/Chapter_10/DCC.stan")
fit_DCC <- sampling(DCC, data = list(N= nrow(x), y=x, Rho_hat=cor(x), K=4, Psi=Psi ,prior=0), chains=1, iter=500)The parameters estimated are similar to Tsay’s (with the exception of the thetas):
print(fit_DCC, names(fit_DCC)[1:15])## Inference for Stan model: DCC.
## 1 chains, each with iter=500; warmup=250; thin=1;
## post-warmup draws per chain=250, total post-warmup draws=250.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## nu 7.68 0.04 0.61 6.68 7.24 7.61 8.06 8.94 255 1.00
## alpha0[1] 0.05 0.00 0.04 0.01 0.02 0.04 0.06 0.17 87 1.00
## alpha0[2] 0.01 0.00 0.00 0.00 0.01 0.01 0.01 0.02 209 1.00
## alpha0[3] 0.01 0.00 0.01 0.00 0.01 0.01 0.01 0.03 138 1.00
## alpha0[4] 0.03 0.00 0.01 0.01 0.02 0.03 0.04 0.06 169 1.01
## alpha1[1] 0.81 0.01 0.13 0.45 0.77 0.84 0.90 0.94 94 1.00
## alpha1[2] 0.94 0.00 0.02 0.90 0.93 0.94 0.95 0.97 196 1.00
## alpha1[3] 0.94 0.00 0.01 0.92 0.94 0.94 0.95 0.97 84 1.00
## alpha1[4] 0.92 0.00 0.01 0.89 0.91 0.92 0.93 0.95 138 1.01
## alpha2[1] 0.03 0.00 0.01 0.01 0.02 0.03 0.04 0.05 113 1.01
## alpha2[2] 0.02 0.00 0.01 0.01 0.02 0.02 0.03 0.04 278 1.00
## alpha2[3] 0.04 0.00 0.01 0.02 0.03 0.04 0.04 0.06 81 1.00
## alpha2[4] 0.06 0.00 0.01 0.04 0.05 0.06 0.07 0.08 141 1.01
## theta1 0.41 0.01 0.10 0.19 0.33 0.43 0.49 0.58 268 1.00
## theta2 0.22 0.01 0.18 0.00 0.08 0.18 0.34 0.64 238 1.00
##
## Samples were drawn using NUTS(diag_e) at Sat Jul 31 15:58:32 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).Additionally, we can plot the estimated volatilities against the rolling estimations:
library(roll)
x <- read.table("data/d-fxsk9904.txt")*100
x <- as.matrix(x); colnames(x) <- c("USEU", "USJP", "IBM", "Dell")
rolling_var <- roll_sd(x, 69)^2
D <- rstan::extract(fit_DCC, pars="D")[[1]];
sigmas <- apply(D, c(2,3,4), mean)
sigmas <- apply(sigmas, 1, diag)^2 %>% t; colnames(sigmas) <- c("USEU", "USJP", "IBM", "Dell")
par(mfrow=c(4,1), mar=c(2,5,1,4))
plot(sigmas[,1], cex.axis=2, cex.lab=2, main="", xlab="", ylab=colnames(sigmas)[1], type="l")
lines(rolling_var[,1], lty=2)
plot(sigmas[,2], cex.axis=2, cex.lab=2, main="", xlab="", ylab=colnames(sigmas)[2], type="l")
lines(rolling_var[,2], lty=2)
plot(sigmas[,3], cex.axis=2, cex.lab=2, main="", xlab="", ylab=colnames(sigmas)[3], type="l")
lines(rolling_var[,3], lty=2)
plot(sigmas[,4], cex.axis=2, cex.lab=2, main="", xlab="", ylab=colnames(sigmas)[4], type="l")
lines(rolling_var[,4], lty=2)
Here we show, the time plots of time-varying correlations between percentage simple returns of four assets:
rolling_cor <- roll_cor(x, x, 69)
Rho_m <- rstan::extract(fit_DCC, pars="Rho_m")[[1]];
USEU_vs_JPUS <- Rho_m[1,,1,2]
IBM_vs_USEU <- Rho_m[1,,1,3]
IBM_vs_JPUS <- Rho_m[1,,2,3]
Dell_vs_USEU <- Rho_m[1,,1,4]
Dell_vs_JPUS <- Rho_m[1,,2,4]
Dell_vs_IBM <- Rho_m[1,,3,4]
par(mfrow=c(6,1), mar=c(2,5,1,4))
plot(USEU_vs_JPUS, cex.axis=2, cex.lab=2, main="", xlab="", ylab="USEU_vs_JPUS", ylim=c(-0.75,0.75), type="l")
lines(rolling_cor[1,2,], lty=2)
plot(IBM_vs_USEU, cex.axis=2, cex.lab=2, main="", xlab="", ylab="IBM_vs_USEU", ylim=c(-0.75,0.75), type="l")
lines(rolling_cor[1,3,], lty=2)
plot(IBM_vs_JPUS, cex.axis=2, cex.lab=2, main="", xlab="", ylab="IBM_vs_JPUS", ylim=c(-0.75,0.75), type="l")
lines(rolling_cor[2,3,], lty=2)
plot(Dell_vs_USEU, cex.axis=2, cex.lab=2, main="", xlab="", ylab="Dell_vs_USEU", ylim=c(-0.75,0.75), type="l")
lines(rolling_cor[1,4,], lty=2)
plot(Dell_vs_JPUS, cex.axis=2, cex.lab=2, main="", xlab="", ylab="Dell_vs_JPUS", ylim=c(-0.75,0.75), type="l")
lines(rolling_cor[2,4,], lty=2)
plot(Dell_vs_IBM, cex.axis=2, cex.lab=2, main="", xlab="", ylab="Dell_vs_IBM", ylim=c(-0.75,0.75), type="l")
lines(rolling_cor[3,4,], lty=2)
Comments
Post a Comment