scripts/Q6_Multivariate_Models.R
In 3schwartz/SpecialeScrAndFun: Script and functions used for master thesis.
#### Constructing_and_sampling_multivariate_distributions ####
## Understanding the construction of multivariate normal distributions,
## normal variance mixtures, elliptical distributions etc.
## Note: This is for educational purposes only; for several of the tasks
## presented below, functions already exist in R packages.
### 1 Basic 'ingredients' ##
## Building a covariance/correlation matrix
## Covariance matrix
A <- matrix(c( 4, 0,
-1, 1), ncol = 2, byrow = TRUE) # Cholesky factor of the ...
(Sigma <- A %*% t(A)) # ... symmetric, positive definite (covariance) matrix Sigma
## Corresponding correlation matrix
P <- outer(1:2, 1:2, Vectorize(function(r, c)
Sigma[r,c]/(sqrt(Sigma[r,r])*sqrt(Sigma[c,c])))) # construct the corresponding correlation matrix
(P. <- cov2cor(Sigma)) # a more elegant solution, see the source of cov2cor()
stopifnot(all.equal(P., P))
## Another option would be as.matrix(Matrix::nearPD(Sigma, corr = TRUE, maxit = 1000)$mat)
## which works differently, though (by finding a correlation matrix close to the
## given matrix in the Frobenius norm) and thus gives a different answer.
## Decomposing a covariance/correlation matrix
## We frequently need to decompose a covariance matrix Sigma (or correlation
## matrix P) as Sigma = AA^T. To this end there are several possibilities.
## The Cholesky decomposition
A. <- t(chol(Sigma)) # the Cholesky factor (lower triangular with real entries > 0)
stopifnot(all.equal(A. %*% t(A.) , Sigma), # checking decomposition
all.equal(A., A)) # checking uniqueness of the Cholesky decomposition
## Other decompositions of Sigma than A %*% t(A) are possible, too
if(FALSE) {
## Eigendecomposition (or spectral decomposition)
eig <- eigen(Sigma) # eigenvalues and eigenvectors
V <- eig$vectors # matrix of eigenvectors
Lambda <- diag(pmax(eig$values, 0))
stopifnot(all.equal(Sigma, V %*% Lambda %*% t(V))) # for real, symmetric matrices
A.eig <- V %*% sqrt(Lambda) %*% t(V)
Sigma.eig <- A.eig %*% t(A.eig)
stopifnot(all.equal(Sigma.eig, Sigma))
## ... but A.. (non-triangular) and A (triangular) are different
## Singular-value decomposition
sv <- svd(Sigma) # singular values, U, V (left/right singular vectors) such that Sigma = U diag(<singular values>) V^T
A.sv <- sv$u %*% sqrt(diag(pmax(sv$d, 0))) %*% t(sv$v)
Sigma.sv <- A.sv %*% t(A.sv)
stopifnot(all.equal(Sigma.sv, Sigma))
}
### 2 Sampling from the multivariate normal distribution ##
## Sampling Z (iid N(0,1))
n <- 1000 # sample size
d <- 2 # dimension
set.seed(271) # set a seed (for reproducibility)
Z <- matrix(rnorm(n * d), ncol = d) # sample iid N(0,1) random variates
mabs <- max(abs(Z))
lim <- c(-mabs, mabs)
plot(Z, xlim = lim, ylim = lim,
xlab = expression(Z[1]), ylab = expression(Z[2]))
## Change the covariance matrix from identity to Sigma
X. <- t(A %*% t(Z)) # ~ N_d(0, Sigma)
xlim <- c(-16, 16)
ylim <- c(-8, 8)
plot(X., xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]))
## Use a shift to get X ~ N_d(mu, Sigma)
mu <- c(1, 3)
X.norm <- rep(mu, each = n) + X.
plot(X.norm, xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]))
## Think about why the rep() and why the t()'s above!
## See also Hofert (2013, "On Sampling from the Multivariate t Distribution")
## Plot (even smaller correlation (larger absolute correlation))
P. <- matrix(c( 1, -0.95,
-0.95, 1), ncol = d, byrow = TRUE)
Sigma. <- outer(1:d, 1:d, Vectorize(function(r, c)
P.[r,c] * (sqrt(Sigma[r,r])*sqrt(Sigma[c,c])))) # construct the corresponding Sigma
## Note: When manually changing the off-diagonal entry of Sigma make sure that
## |Sigma.[1,2]| <= sqrt(Sigma.[1,1]) * sqrt(Sigma.[2,2]) still holds
A. <- t(chol(Sigma.)) # new Cholesky factor
X.norm. <- rep(mu, each = n) + t(A. %*% t(Z)) # generate the sample
plot(rbind(X.norm, X.norm.), xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("black", "royalblue3"), each = n))
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
legend = c(expression(N(mu,Sigma)),
expression("With larger absolute correlation")))
## Plot (switch sign of correlation)
Sigma. <- Sigma * matrix(c( 1, -1,
-1, 1), ncol = d, byrow = TRUE)
A. <- t(chol(Sigma.))
X.norm. <- rep(mu, each = n) + t(A. %*% t(Z))
plot(rbind(X.norm, X.norm.), xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("black", "royalblue3"), each = n))
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
legend = c(expression(N(mu,Sigma)),
expression("Correlation sign switched")))
## Note: This has all been done based on the same Z (and we also recycle it below)!
### 3 Sampling from the multiv. t distribution (and other normal variance mixtures)
## We recycle the sample with positive correlation
X.norm <- X.norm.
A <- A.
Sigma <- A %*% t(A)
### 3.1 A sample from a multivariate t_nu distribution ##
nu <- 3 # degrees of freedom
set.seed(314)
W <- 1/rgamma(n, shape = nu/2, rate = nu/2) # mixture variable for a multiv. t_nu distribution
plot(W, log = "y")
X.t <- rep(mu, each = n) + sqrt(W) * t(A %*% t(Z))
xlim <- c(-40, 80)
ylim <- c(-20, 60)
plot(rbind(X.t, X.norm), xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("maroon3", "black"), each = n))
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "maroon3"),
legend = c(expression(N(mu,Sigma)),
expression(italic(t)[nu](mu,Sigma))))
## For more about sampling from the multivariate t distribution (especially
## what can go wrong!), see Hofert (2013, "On Sampling from the Multivariate t Distribution")
### 3.2 A 2-point distribution for W (= two 'overlaid' multiv. normals) ##
## Let's change the distribution of W to a 2-point distribution with
## W = 1 with prob. 0.5 and W = 2 with prob. 0.5. What we will see are two 'overlaid'
## normal distributions (as each of W = 1 and W = 2 corresponds to one).
W.binom <- 1 + rbinom(n, size = 1, prob = 0.5)
cols <- rep("royalblue3", n)
cols[W.binom == 2] <- "maroon3"
X.W.binom <- rep(mu, each = n) + sqrt(W.binom) * t(A %*% t(Z))
plot(rbind(X.t, X.W.binom), xlab = expression(X[1]), ylab = expression(X[2]),
col = c(rep("black", n), cols))
legend("topright", bty = "n", pch = rep(1, 3), col = c("royalblue3", "maroon3", "black"),
legend = c(expression(N(mu,W*Sigma)~"for W = 1"),
expression(N(mu,W*Sigma)~"for W = 2"),
expression(italic(t)[nu](mu,Sigma))))
## => With probability 0.5 we sample from a normal variance mixture with W = 1
## and with probability 0.5 we sample from one with W = 2. By using a df for
## W with infinite upper endpoint, we can reach further out in the tails than
## with any multivariate normal distribution by overlaying normals
## with different (unbounded) covariance matrices (and this is what is
## creating heavier tails for the t distribution than the normal).
### 3.3 An *uncorrelated* multivariate t_nu distribution vs an independent one
A. <- diag(sqrt(diag(Sigma))) # Cholesky factor of the diagonal matrix containing the variances
X.t.uncor <- rep(mu, each = n) + sqrt(W) * t(A. %*% t(Z)) # uncorrelated t samples
X.t.ind <- rep(mu, each = n) + matrix(rt(n * d, df = nu), ncol = d) # independent t samples
plot(rbind(X.t.uncor, # uncorrelated t; dependence only introduced by sqrt(W) => shifted spherical distribution
X.t, # full t; dependence introduced by sqrt(W) *and* correlation
X.t.ind), # independence
xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("royalblue3", "maroon3", "black"), each = n))
legend("topright", bty = "n", pch = rep(1, 3), col = c("black", "royalblue3", "maroon3"),
legend = c(expression("Independent"~italic(t)[nu]),
expression("Uncorrelated"~italic(t)[nu]~"(shifted spherical)"),
expression(italic(t)[nu](mu,Sigma))))
## => The uncorrelated (but not independent) samples show more mass in the joint
## tails than the independent samples. The (neither uncorrelated nor independent)
## t samples show the most mass in the joint tails.
### 3.4 A normal *mean*-variance mixture ##
## Now let's also 'mix' the mean (replace mu by mu(W)), to get a normal
## mean-variance mixture; here: choosing between two different locations mu
## depending on W.
mu. <- t(sapply(W.binom, function(w) mu + if(w == 1) 0 else c(15, 0)))
X.mean.var <- mu. + sqrt(W.binom) * t(A %*% t(Z))
plot(rbind(X.W.binom, X.mean.var), xlab = expression(X[1]), ylab = expression(X[2]),
col = c(rep("black", n), cols))
legend("topleft", bty = "n", pch = rep(1, 3), col = c("black", "royalblue3", "maroon3"),
legend = c(expression(N(mu,W*Sigma)~"for W"%in%"{1,2}"),
"Mean-var. mix for W = 1",
"Mean-var. mix for W = 2"))
## => Not only do the red points show more variation, they *also* have a different
## location; clearly, this sample (blue + red points together) is not
## elliptically distributed anymore.
### 4 Sampling elliptical distributions ##
### 4.1 Sample from a multivariate t_nu via elliptical distributions ##
## Sample from the uniform distribution on the unit sphere (spherical distribution)
set.seed(271)
Z <- matrix(rnorm(n * d), ncol = d)
S <- Z/sqrt(rowSums(Z^2)) # Z/||Z||
lim <- c(-1.3, 1.3)
par(pty = "s")
plot(S, xlim = lim, ylim = lim,
xlab = expression(S[1]), ylab = expression(S[2]))
## Add scale (elliptical distribution) by computing X = RAS (= ARS)
rho <- 0.8 # correlation
P <- matrix(c( 1, rho,
rho, 1), ncol = 2) # correlation matrix
A <- chol(P) # Cholesky factor (here: upper triangular matrix for efficiency reasons)
X. <- t(A %*% t(S)) # (n, d) matrix
par(pty = "s")
plot(X., xlim = lim, ylim = lim,
xlab = expression(X[1]), ylab = expression(X[2]))
## Add radial part of a t_nu distribution (spherical distribution)
nu <- 5
R <- sqrt(d * rf(n, df1 = d, df2 = nu))
X.. <- R * X.
ran <- c(min(X.., -max(X..)), max(X.., -min(X..)))
par(pty = "s")
plot(X.., xlim = ran, ylim = ran, xlab = expression(X[1]), ylab = expression(X[2]))
## Add location (elliptical distribution) by computing X = mu + RAS
mu <- c(ran[2]-max(X..[,1]), ran[2]-max(X..[,2])) # push sample towards the upper right corner
X <- rep(mu, each = n) + X..
par(pty = "s")
plot(X, xlim = ran, ylim = ran, xlab = expression(X[1]), ylab = expression(X[2]))
### 4.2 A 2-point distribution for R ##
R <- 1 + rbinom(n, size = 1, prob = 2/3) # P(R = 1) = 1/3, P(R = 2) = 2/3
cols <- rep("royalblue3", n)
cols[R == 2] <- "maroon3" # prob. 2/3
X <- rep(mu, each = n) + t(A %*% t(R*S)) # same mu, A, S
par(pty = "s")
plot(X, xlab = expression(X[1]), ylab = expression(X[2]), col = cols)
legend("bottomright", bty = "n", pch = c(1,1), col = c("royalblue3", "maroon3"),
legend = c("Elliptical for R = 1", "Elliptical for R = 2"))
## ... in other words, the radial part R scales each point on the ellipse AS
## to lie on another ellipse. For continuously distributed R, these ellipses
## are all overlaid (and thus not visible anymore), but for discrete R, we
## may see them (here: two ellipses).
#### Fitting_a_multivariate_normal_and_t ####
## Fitting a multivariate normal and t distribution to negative log-returns of
## the five components of the Dow Jones index
library(xts) # for time series manipulation
library(QRM) # for fit.mst()
library(mvtnorm) # for sampling from a multivariate normal or t distribution
library(qrmdata) # for Dow Jones constituents data
library(qrmtools) # for returns()
## Load and extract the data we work with and plot
data(DJ_const)
str(DJ_const)
S <- DJ_const['2000-01-01/', c("AAPL", "BA", "INTC", "IBM", "NKE")]
## => We work with Apple, Boeing, Intel, IBM, Nike since 2000-01-01 here
plot.zoo(S, xlab = "Time t")
## Build and plot -log-returns
X <- -returns(S) # compute -log-returns
pairs(as.matrix(X), main = "Scatter plot matrix of risk-factor changes", gap = 0, pch = ".")
plot.zoo(X, xlab = "Time t", main = "")
## For illustration purposes, we treat this data as realizations of iid
## five-dimensional random vectors.
## Fitting a multivariate normal distribution to X and simulating from it
mu <- colMeans(X) # estimated location vector
Sigma <- cov(X) # estimated scale matrix
stopifnot(all.equal(Sigma, var(X)))
P <- cor(X) # estimated correlation matrix
stopifnot(all.equal(P, cov2cor(Sigma)))
n <- nrow(X) # sample size
set.seed(271)
X.norm <- rmvnorm(n, mean = mu, sigma = Sigma) # N(mu, Sigma) samples
## Fitting a multivariate t distribution to X
fit <- fit.mst(X, method = "BFGS") # fit a multivariate t distribution
X.t <- rmvt(n, sigma = as.matrix(fit$Sigma), df = fit$df, delta = fit$mu) # t_nu samples
## Plot (sample from fitted t (red), original sample (black), sample from fitted normal (blue))
dat <- rbind(t = X.t, original = as.matrix(X), norm = X.norm)
cols <- rep(c("maroon3", "black", "royalblue3"), each = n)
pairs(dat, gap = 0, pch = ".", col = cols)
## Pick out one pair (to better see that multivariate t fits better)
plot(dat[,1:2], col = cols) # => the multivariate normal generates too few extreme losses!
legend("bottomright", bty = "n", pch = rep(1, 3), col = c("black", "royalblue3", "maroon3"),
legend = c("-Log-returns", expression("Simulated fitted"~N(mu,Sigma)),
expression("Simulated fitted"~italic(t)[nu](mu,Sigma))))
#### Multivariate_normal ####
## Libraries required
library(QRM)
library(xts)
library(mvtnorm)
library(qrmdata)
library(copula)
## The bivariate normal distribution
BiDensPlot(func = dmnorm, xpt = c(-4, 4), ypts = c(-4, 4),
mu = c(0, 0), Sigma = equicorr(2, -0.7))
## Some trivariate simulated data
rho <- 0.7
P <- matrix(rho,ncol = 3,nrow = 3)+(1-rho)*diag(3)
P
data <- rmvnorm(1000, sigma = P)
## ML estimation
n <- dim(data)[1]
mu.hat <- colMeans(data)
Sigma.hat <- var(data)*(n-1)/n
P.hat <- cor(data)
## Note that var(data) gives unbiased estimator which is not MLE
ll.max <- sum(dmvnorm(data, mean = mu.hat, sigma = Sigma.hat, log = TRUE))
ll.max
## Some 10-d simulated data
rho = 0.6
P <- matrix(rho, ncol = 10, nrow = 10)+(1-rho)*diag(10)
P
data <- rmvnorm(1000, sigma = P)
## Test of marginal normality
apply(data, 2, shapiro.test)
## Tests of joint normality
MardiaTest(data)
jointnormalTest(data)
## Dow Jones data
data("DJ_const")
Sdata <- DJ_const['2000-01-01/',1:10]
Xdata <- diff(log(Sdata))[-1,]
plot.zoo(Xdata)
Xdata.w <- apply.weekly(Xdata,FUN = colSums)
Xdata.m <- apply.monthly(Xdata,FUN = colSums)
Xdata.q <- apply.quarterly(Xdata,FUN = colSums)
## Test of marginal normality
apply(Xdata, 2, shapiro.test)
apply(Xdata.w, 2, shapiro.test)
apply(Xdata.m, 2, shapiro.test)
apply(Xdata.q, 2, shapiro.test)
## Test for joint normality
MardiaTest(as.matrix(Xdata))
jointnormalTest(as.matrix(Xdata))
MardiaTest(as.matrix(Xdata.w))
jointnormalTest(as.matrix(Xdata.w))
MardiaTest(as.matrix(Xdata.m))
jointnormalTest(as.matrix(Xdata.m))
MardiaTest(as.matrix(Xdata.q))
jointnormalTest(as.matrix(Xdata.q))
## Constructing data which are only marginally normal
set.seed(133)
G.cop <- archmCopula("Gumbel", param = 3, dim = 3)
U.G <- rCopula(1000, copula = G.cop)
data = qnorm(U.G)
pairs(data)
apply(data,2,shapiro.test)
MardiaTest(data)
jointnormalTest(data)
#### Multivariate_Student ####
library(xts)
library(mvtnorm)
library(QRM) # for fit.mst()
library(qrmdata)
## Simulate multivariate Student
rho <- 0.7
d <- 3
P <- matrix(rho, ncol = d, nrow = d) + (1-rho)*diag(d)
nu <- 10
data <- rmvt(2000, sigma = P, df = nu)
pairs(data)
## Test of marginal normality
apply(data, 2, shapiro.test)
## Tests of joint normality
MardiaTest(data)
jointnormalTest(data)
## QQplots against Student t
qplotfunc <- function(data){
qqplot(qt(ppoints(data), df = nu), data, xlab = "Theoretical", ylab = "Empirical")
qqline(data, dist = function(p) qt(p, df = nu))
}
par(mfrow = c(1,3))
apply(data, 2, qplotfunc)
par(mfrow = c(1, 1))
## QQplot of Mahalanobis distances (scaled by dimension d)
## against an F(d,nu) distribution
Ddata <- mahalanobis(data, center = FALSE, cov = var(data))/d
qqplot(qf(ppoints(Ddata), df1 = d, df2 = nu), Ddata,
xlab = "Theoretical", ylab = "Empirical")
qqline(Ddata, dist = function(p) qf(p,df1 = d, df2 = nu))
## Estimate multivariate Student
## Dow Jones Data
data("DJ_const")
Sdata <- DJ_const['2000-01-01/',1:10]
Xdata <- diff(log(Sdata))[-1,]
Xdata.w <- apply.weekly(Xdata, FUN = colSums)
Xdata.m <- apply.monthly(Xdata, FUN = colSums)
## Write function for fitting multivariate normal by MLE
fit.mvnorm <- function(data){
n <- dim(data)[1]
mu.hat <- colMeans(data)
Sigma.hat <- var(data)*(n-1)/n
ll.max <- sum(dmvnorm(data, mean = mu.hat, sigma = Sigma.hat, log = TRUE))
list(mu = mu.hat, Sigma = Sigma.hat, ll.max = ll.max)
}
## Fit normal
mod.w.norm = fit.mvnorm(as.matrix(Xdata.w))
mod.w.norm
## Fit Student t
mod.w.t = fit.mst(as.matrix(Xdata.w), method = "BFGS")
mod.w.t
names(mod.w.t)
mod.w.t$df
## Compare
mod.w.t$ll.max - mod.w.norm$ll.max
## Fit normal
mod.m.norm = fit.mvnorm(as.matrix(Xdata.m))
mod.m.norm
## Fit Student t
mod.m.t = fit.mst(as.matrix(Xdata.m), method = "BFGS")
mod.m.t
## Compare
mod.m.t$ll.max - mod.m.norm$ll.max
#### Multivariate_GH ####
## Libraries required
library(xts)
library(qrmdata)
library(ghyp)
## Dow Jones Data
data("DJ_const")
Sdata <- DJ_const['2000-01-01/',1:10]
Xdata <- diff(log(Sdata))[-1,]
Xdata.w <- apply.weekly(Xdata, FUN = colSums) # weekly data
plot.zoo(Xdata.w)
## Fit elliptical models
mod.NIG.sym <- fit.NIGmv(Xdata.w, symmetric = TRUE, nit = 10000)
summary(mod.NIG.sym)
mod.T.sym <- fit.tmv(Xdata.w, symmetric = TRUE, nit = 10000)
summary(mod.T.sym)
mod.GHYP.sym <- fit.ghypmv(Xdata.w, symmetric = TRUE, nit = 10000)
summary(mod.GHYP.sym)
## Fit skewed models
mod.NIG.skew <- fit.NIGmv(Xdata.w, symmetric = FALSE, nit = 10000)
summary(mod.NIG.skew)
mod.T.skew <- fit.tmv(Xdata.w, symmetric = FALSE, nit = 10000)
summary(mod.T.skew)
mod.GHYP.skew <- fit.ghypmv(Xdata.w, symmetric = FALSE, nit = 10000)
summary(mod.GHYP.skew)
mod.NIG.sym@llh
mod.NIG.skew@llh
mod.T.sym@llh
mod.T.skew@llh
mod.GHYP.sym@llh
mod.GHYP.skew@llh
## Hypothesis: elliptical T is good enough
LRstat <- 2*(mod.T.skew@llh-mod.T.sym@llh)
1-pchisq(LRstat,10) # => H0 not rejected
## Hypothesis: elliptical ghyp is good enough
LRstat <- 2*(mod.GHYP.skew@llh-mod.GHYP.sym@llh)
1-pchisq(LRstat,10) # => H0 not rejected
## Hypothesis: skewed t no worse than skewed ghyp
LRstat <- 2*(mod.GHYP.skew@llh-mod.T.skew@llh)
1-pchisq(LRstat,1) # => H0 rejected at 5% level
## Hypothesis: elliptical t no worse than elliptical ghyp
LRstat <- 2*(mod.GHYP.sym@llh-mod.T.sym@llh)
1-pchisq(LRstat,1) # => H0 rejected at 5% level
## Find skewness parameters
getSlots("mle.ghyp")
mod.GHYP.skew@gamma
## Make picture of two returns
Xdata <- Xdata.w[,c(5,8)]
plot(as.numeric(Xdata[,1]), as.numeric(Xdata[,2]),
xlab = names(Xdata)[1], ylab = names(Xdata)[2])
mod <- fit.ghypmv(Xdata, symmetric = FALSE, nit = 10000)
summary(mod)
xvals <- seq(from = min(Xdata[,1]), to = max(Xdata[,2]), length.out = 50)
yvals <- seq(from = min(Xdata[,2]), to = max(Xdata[,2]), length.out = 50)
outer.func <- function(x,y) dghyp(cbind(x,y), mod)
zvals <- outer(xvals, yvals, outer.func)
zvals.un <- unique(as.numeric(zvals))
contour(xvals, yvals, zvals, add = TRUE, col = 2,
levels = seq(from = quantile(zvals.un,0.6), to = quantile(zvals.un,0.9), length = 10))
#### Fundamental_factor_model ####
## Fitting a fundamental factor model using industry sector information
library(xts) # for time series manipulation
library(qrmdata) # for SP500 (constituents) data
library(qrmtools) # for plot_NA() and plot_matrix()
### 1 Data preparation ##
## Load and extract the data we work with (all available since 1995)
## Index
data(SP500_const) # index data
SP500.const <- SP500_const['1995-01-01/',] # all since 1995
dim(SP500.const)
NA_plot(SP500.const[,1:100]) # => many stocks have missing values (look at first 100)
SP500.const <- SP500.const[, colSums(is.na(SP500.const)) <= 0.1 * nrow(SP500.const)] # omit columns with more than 10% NA
SP500.const <- na.fill(SP500.const, fill = "extend") # fill the remaining NAs
dim(SP500.const)
## Caution: only execute next block if you want to see all the data plotted!!!
par(ask = TRUE)
nplots <- ceiling(dim(SP500.const)[2]/20)
for (i in 1:nplots){
r1 <- 20*(i-1)+1
r2 <- min(20*i,dim(SP500.const)[2])
plot.zoo(SP500.const[,r1:r2], xlab = "Time t", main = "SP500 Constituents")
}
par(ask = FALSE)
## Remove certain stocks
## AIG, C, ETFC (badly affected by financial crisis)
## T, CTL, FTR, VZ (only 4 telecommunications companies so this sector is thin)
## CMCSK (moves together with CMCSA)
## FOX (moves together with FOXA)
ignore <- list("AIG","T","CTL","C","CMCSK","ETFC","FTR","FOX","VZ")
which(names(SP500.const) %in% ignore)
SP500.const <- SP500.const[,-which(names(SP500.const) %in% ignore)]
dim(SP500.const)
## Build and plot log-returns
X.const <- diff(log(SP500.const))[-1,] # compute -log-returns
## We use monthly data here as basis (and compute monthly log-returns)
X <- apply.monthly(X.const, FUN = colSums) # 252 by 367 matrix
plot.zoo(X[,1:20], type = "h", xlab = "Time t",
main = "Monthly risk-factor changes (log-returns) of SP500 constituents")
summary(SP500_const_info)
firms <- as.character(SP500_const_info$Ticker)
## Make two firm names consistent
firms[firms == "BRK-B"] <- "BRK.B"
firms[firms == "BF-B"] <- "BF.B"
sectors <- as.character(SP500_const_info$Sector)
## Construct a factor X.sectors and tabulate levels
X.sectors <- sectors[which(firms %in% colnames(X))]
X.sectors <- as.factor(X.sectors)
table(X.sectors)
levels(X.sectors) <- c("Con-Disc.","Con-Stap.","Energy",
"Financials","Health","Industrials","IT",
"Materials","Utilities")
### 2 Model fitting ##
## Fit a cross-sectional regression model X_t = B*F_t + eps
## A regression is fitted at each time point
## No intercept is required (-1)
mod <- lm(t(X) ~ X.sectors -1)
B <- model.matrix(mod)
## Factors can be constructed from coefficient matrix
coef.mat <- t(coef(mod))
dimnames(coef.mat)[[2]] <- levels(X.sectors)
F <- xts(coef.mat,time(X))
plot.zoo(F,type = "h")
## Calculate error variances for each stock
eps <- t(residuals(mod))
eps.variances <- diag(var(eps))
## Display them from smallest to largest
sort(eps.variances)
## Use generalized least squares to obtain better estimates
## See documentation of lm for more details
mod2 <- lm(t(X) ~ X.sectors -1, weights = 1/eps.variances)
coef.mat <- t(coef(mod2))
dimnames(coef.mat)[[2]] <- levels(X.sectors)
F <- xts(coef.mat,time(X))
plot.zoo(F,type = "h")
cor(F)
eps <- t(residuals(mod))
## This matrix should be near-diagonal for perfect factor model
## Of course this won't be the case in practice
cor.eps <- cor(eps)
## Is Cor(eps) (roughly) diagonal? Check upper corner
diag(cor.eps) <- NA
matrix_plot(cor.eps[1:75,1:75])
## Are the errors uncorrelated with the factors?
cor.eps.F <- cor(eps, F)
summary(cor.eps.F) # => there are some large correlations
## Construct the implied covariance and correlation matrix
Ups <- cov(eps) # Upsilon (covariance matrix of epsilon)
Omega <- as.matrix(cov(F)) # Omega (covariance matrix of F; systematic risk)
Sigma <- B %*% Omega %*% t(B) + diag(diag(Ups)) # Cov(X)
P <- cov2cor(Sigma) # Cor(X)
## Look at discrepancies between the factor model correlation matrix and the
## sample correlation matrix; check upper corner
err <- P-cor(X)
matrix_plot(err[1:75,1:75])
#### Macroeconomic_factor_model ####
## Fitting a factor model
library(xts) # for time series manipulation
library(qrmdata) # for Dow Jones (constituents) data
library(qrmtools) # for plot_NA() and plot_matrix()
### 1 Data preparation ##
## Load and extract the data we work with (all available since 1990) and plot
## Index
data(DJ) # index data
DJ. <- DJ['1990-01-01/'] # all since 1990
plot.zoo(DJ., xlab = "Time t", ylab = "Dow Jones Index")
## Constituents
data(DJ_const) # constituents data
DJ.const <- DJ_const['1990-01-01/',] # all since 1990
NA_plot(DJ.const) # => use all but the two columns with lots of NAs
DJ.const <- DJ.const[, colSums(is.na(DJ.const)) <= 0.1 * nrow(DJ.const)] # omit columns with more than 10% NA
DJ.const <- na.fill(DJ.const, fill = "extend") # fill the remaining NAs
plot.zoo(DJ.const, xlab = "Time t", main = "Dow Jones Constituents")
## Build and plot log-returns
## Index
X. <- diff(log(DJ.))[-1,] # compute -log-returns
plot.zoo(X., xlab = "Time t", ylab = expression(X[t]), main = "Risk-factor changes (log-returns) of Dow Jones index")
## Constituents
X.const <- diff(log(DJ.const))[-1,] # compute -log-returns
if(FALSE) # more time-consuming
pairs(as.matrix(X.const), gap = 0, pch = ".",
main = "Scatter plot matrix of risk-factor changes (log-returns) of Dow Jones constituents")
plot.zoo(X.const, xlab = "Time t", main = "Risk-factor changes (log-returns) of Dow Jones constituents")
## We use monthly data here as basis (and compute monthly log-returns)
X <- apply.monthly(X.const, FUN = colSums) # (312, 28)-matrix
plot.zoo(X, type = "h", xlab = "Time t", main = "Monthly risk-factor changes (log-returns) of Dow Jones constituents")
F <- apply.monthly(X., FUN = colSums)
plot(F, type = "h", xlab = "Time t", ylab = expression(X[t]), main = "Monthly risk-factor changes (log-returns) of Dow Jones index")
### 2 Model fitting ##
## Fit a multivariate regression model X = a + B*F + eps
(res <- lm(X ~ F)) # more details via summary()
## Get parameter estimates
names(res)
par.ests <- coefficients(res)
a <- par.ests[1,]
B <- par.ests[2,]
# tail(X[,1:3])
# tail(predict(res)[,1:3])
# forecast::forecast(res, newdata = F[length(F)], h = 1)
## Get residuals and estimate their correlation matrix
eps <- resid(res) # (312, 28)-matrix
cor.eps <- cor(eps)
## Is Cor(eps) (roughly) diagonal?
matrix_plot(cor.eps, at = seq(-1, 1, length.out = 200)) # => yes (as required)
## Are the errors uncorrelated with the factors?
cor.eps.F <- cor(eps, F) # 28 correlations (idiosyncratic risk)
summary(cor.eps.F) # => yes (as required)
## Construct the implied covariance and correlation matrix
Ups <- cov(eps) # Upsilon (covariance matrix of epsilon)
Omega <- as.matrix(cov(F)) # Omega (covariance matrix of F; systematic risk)
Sigma <- B %*% Omega %*% t(B) + diag(diag(Ups)) # Cov(X); (28, 28)-matrix
rownames(Sigma) <- colnames(Sigma)
P <- cov2cor(Sigma) # Cor(X)
## Look at discrepancies between the factor model correlation matrix and the
## sample correlation matrix
err <- P-cor(X)
matrix_plot(err, at = seq(-1, 1, length.out = 200))
#### PCA_factor_model ####
## Fitting a PCA factor model
library(xts) # for time series manipulation
library(qrmdata) # for Dow Jones (constituents) data
library(qrmtools) # for plot_NA() and plot_matrix()
### 1 Data preparation ###
## Load and extract the data we work with (all available since 1990) and plot
## Index
data(DJ) # index data
DJ. <- DJ['1990-01-01/'] # all since 1990
plot.zoo(DJ., xlab = "Time t", ylab = "Dow Jones Index")
## Constituents
data(DJ_const) # constituents data
DJ.const <- DJ_const['1990-01-01/',] # all since 1990
NA_plot(DJ.const) # => use all but the two columns with lots of NAs
DJ.const <- DJ.const[, colSums(is.na(DJ.const)) <= 0.1 * nrow(DJ.const)] # omit columns with more than 10% NA
DJ.const <- na.fill(DJ.const, fill = "extend") # fill the remaining NAs
plot.zoo(DJ.const, xlab = "Time t", main = "Dow Jones Constituents")
## Build and plot log-returns
## Index
X. <- diff(log(DJ.))[-1,] # compute -log-returns
plot.zoo(X., xlab = "Time t", ylab = expression(X[t]),
main = "Risk-factor changes (log-returns) of Dow Jones index")
## Constituents
X.const <- diff(log(DJ.const))[-1,] # compute -log-returns
if(FALSE) # more time-consuming
pairs(as.matrix(X.const), gap = 0, pch = ".",
main = "Scatter plot matrix of risk-factor changes (log-returns) of Dow Jones constituents")
plot.zoo(X.const, xlab = "Time t", main = "Risk-factor changes (log-returns) of Dow Jones constituents")
## We use monthly data here as basis (and compute monthly log-returns)
X <- apply.monthly(X.const, FUN = colSums) # (312, 28)-matrix
plot.zoo(X, type = "h", xlab = "Time t",
main = "Monthly risk-factor changes (log-returns) of Dow Jones constituents")
Xindex <- apply.monthly(X., FUN = colSums)
plot(Xindex, type = "h", xlab = "Time t", ylab = expression(X[t]),
main = "Monthly risk-factor changes (log-returns) of Dow Jones index")
### 2 Model fitting ###
## Principal component analysis
PCA <- princomp(X) # determine the principal components
PCA. <- prcomp(X) # another option
summary(PCA)
plot(PCA) # show important components
loadings(PCA) # factor loadings (representing how much a factor explains a variable)
## Working with the principal components
mu <- PCA$center # estimated centers
Y <- PCA$scores # estimated principal components of X
nprin <- 3 # number of important principal components
Y1 <- xts(Y[,1:nprin], time(X)) # grab out the first so-many important principal components of X
plot.zoo(Y1, type = "h", xlab = "Time", ylab = paste("Component", 1:nprin),
main = "Principal components of X")
cor(Y1)
Y2 <- Y[,(nprin+1):ncol(Y)] # ignored principal components of X
## Interpreting first principal component
plot.zoo(merge(Xindex,Y1[,1]),main=" ") # compare DJ index returns and first PC
plot(as.numeric(Xindex),-as.numeric(Y1[,1]),xlab="Index",ylab="PC1 (neg)")
## Reconstruct X from Y1 and Y2
G <- unclass(loadings(PCA)) # compute G, see McNeil, Frey, Embrechts (2015, (6.64))
G1 <- G[,1:nprin] # compute G_1, see McNeil, Frey, Embrechts (2015, (6.65))
G2 <- G[,(nprin+1):ncol(G)]
eps <- G2 %*% t(Y2) # compute epsilon
eps.cor <- cor(t(eps)) # compute correlations of epsilons
diag(eps.cor) <- NA # neglect diagonal elements
matrix_plot(eps.cor) # not perfectly diagonal
X. <- t(mu + G1 %*% t(Y1) + eps)
err <- X.-X
head(err) # check that differences are zero vectors
summary(as.vector(err))
## Note: We could fit a univariate GARCH model to each PCA factor.
## This is called PC-GARCH model.
#### Projections_of_uniform_on_sphere ####
## Orthogonal projections of the uniform distribution on the unit sphere S_{d-1}
## in d dimensions down by two dimensions are uniform in the (d-2)-ball;
## see http://www.math.univ-toulouse.fr/~letac/Archimede.pdf
## Let's check if that's the case!
library(rgl)
## Sample from the uniform distribution on the unit d-sphere S_{d-1}
n <- 5e3 # sample size
d <- 4 # maximal dimension investigated
set.seed(271) # for reproducibility
Z <- matrix(rnorm(n*d), ncol = d) # generate independent N(0,1)
### 1 Projections from d = 3 to d = 2 ##
S2 <- Z[,1:3]/sqrt(rowSums(Z[,1:3]^2)) # ~ U(S_2)
pairs(S2, gap = 0, pch = ".") # => does not look uniform (fine)
### 2 Projections from d = 3 to d = 1 ##
plot(S2[,1], ylab = expression(S[1]), cex = 0.5) # => looks uniform in [0,1]
plot(S2[,2], ylab = expression(S[2]), cex = 0.5) # => looks uniform in [0,1]
plot(S2[,3], ylab = expression(S[3]), cex = 0.5) # => looks uniform in [0,1]
### 3 Projections from d = 4 to d = 3 ##
S3 <- Z[,1:4]/sqrt(rowSums(Z[,1:4]^2)) # ~ U(S_3)
colnames(S3) <- paste0("S", 1:4)
plot3d(S3[,c(2,3,4)]) # => looks uniform in the ball in 3d but difficult to judge
S3. <- S3[0.2 < S3[,2] & S3[,2] <= 0.3, c(2,3,4)] # pick out a slice (0.2, 0.3]
plot(S3.) # slice looks uniform
plot3d(S3[,c(1,3,4)])
plot3d(S3[,c(1,2,4)])
plot3d(S3[,c(1,2,3)])
## Note: Although not covered by the above result, it *seems* that projections
## from U(S_3) to IR^3 are also uniform in the 3-ball. More investigations
## necessary.
### 4 Projections from d = 4 to d = 2 ##########################################
pairs(S3, gap = 0, pch = ".", # => looks (pairwise) uniform in the unit disk
labels = as.expression(sapply(1:4, function(j) bquote(S[.(j)]))))
3schwartz/SpecialeScrAndFun documentation built on May 4, 2019, 6:29 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#### Constructing_and_sampling_multivariate_distributions ####
## Understanding the construction of multivariate normal distributions,
## normal variance mixtures, elliptical distributions etc.
## Note: This is for educational purposes only; for several of the tasks
## presented below, functions already exist in R packages.
### 1 Basic 'ingredients' ##
## Building a covariance/correlation matrix
## Covariance matrix
A <- matrix(c( 4, 0,
-1, 1), ncol = 2, byrow = TRUE) # Cholesky factor of the ...
(Sigma <- A %*% t(A)) # ... symmetric, positive definite (covariance) matrix Sigma
## Corresponding correlation matrix
P <- outer(1:2, 1:2, Vectorize(function(r, c)
Sigma[r,c]/(sqrt(Sigma[r,r])*sqrt(Sigma[c,c])))) # construct the corresponding correlation matrix
(P. <- cov2cor(Sigma)) # a more elegant solution, see the source of cov2cor()
stopifnot(all.equal(P., P))
## Another option would be as.matrix(Matrix::nearPD(Sigma, corr = TRUE, maxit = 1000)$mat)
## which works differently, though (by finding a correlation matrix close to the
## given matrix in the Frobenius norm) and thus gives a different answer.
## Decomposing a covariance/correlation matrix
## We frequently need to decompose a covariance matrix Sigma (or correlation
## matrix P) as Sigma = AA^T. To this end there are several possibilities.
## The Cholesky decomposition
A. <- t(chol(Sigma)) # the Cholesky factor (lower triangular with real entries > 0)
stopifnot(all.equal(A. %*% t(A.) , Sigma), # checking decomposition
all.equal(A., A)) # checking uniqueness of the Cholesky decomposition
## Other decompositions of Sigma than A %*% t(A) are possible, too
if(FALSE) {
## Eigendecomposition (or spectral decomposition)
eig <- eigen(Sigma) # eigenvalues and eigenvectors
V <- eig$vectors # matrix of eigenvectors
Lambda <- diag(pmax(eig$values, 0))
stopifnot(all.equal(Sigma, V %*% Lambda %*% t(V))) # for real, symmetric matrices
A.eig <- V %*% sqrt(Lambda) %*% t(V)
Sigma.eig <- A.eig %*% t(A.eig)
stopifnot(all.equal(Sigma.eig, Sigma))
## ... but A.. (non-triangular) and A (triangular) are different
## Singular-value decomposition
sv <- svd(Sigma) # singular values, U, V (left/right singular vectors) such that Sigma = U diag(<singular values>) V^T
A.sv <- sv$u %*% sqrt(diag(pmax(sv$d, 0))) %*% t(sv$v)
Sigma.sv <- A.sv %*% t(A.sv)
stopifnot(all.equal(Sigma.sv, Sigma))
}
### 2 Sampling from the multivariate normal distribution ##
## Sampling Z (iid N(0,1))
n <- 1000 # sample size
d <- 2 # dimension
set.seed(271) # set a seed (for reproducibility)
Z <- matrix(rnorm(n * d), ncol = d) # sample iid N(0,1) random variates
mabs <- max(abs(Z))
lim <- c(-mabs, mabs)
plot(Z, xlim = lim, ylim = lim,
xlab = expression(Z[1]), ylab = expression(Z[2]))
## Change the covariance matrix from identity to Sigma
X. <- t(A %*% t(Z)) # ~ N_d(0, Sigma)
xlim <- c(-16, 16)
ylim <- c(-8, 8)
plot(X., xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]))
## Use a shift to get X ~ N_d(mu, Sigma)
mu <- c(1, 3)
X.norm <- rep(mu, each = n) + X.
plot(X.norm, xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]))
## Think about why the rep() and why the t()'s above!
## See also Hofert (2013, "On Sampling from the Multivariate t Distribution")
## Plot (even smaller correlation (larger absolute correlation))
P. <- matrix(c( 1, -0.95,
-0.95, 1), ncol = d, byrow = TRUE)
Sigma. <- outer(1:d, 1:d, Vectorize(function(r, c)
P.[r,c] * (sqrt(Sigma[r,r])*sqrt(Sigma[c,c])))) # construct the corresponding Sigma
## Note: When manually changing the off-diagonal entry of Sigma make sure that
## |Sigma.[1,2]| <= sqrt(Sigma.[1,1]) * sqrt(Sigma.[2,2]) still holds
A. <- t(chol(Sigma.)) # new Cholesky factor
X.norm. <- rep(mu, each = n) + t(A. %*% t(Z)) # generate the sample
plot(rbind(X.norm, X.norm.), xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("black", "royalblue3"), each = n))
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
legend = c(expression(N(mu,Sigma)),
expression("With larger absolute correlation")))
## Plot (switch sign of correlation)
Sigma. <- Sigma * matrix(c( 1, -1,
-1, 1), ncol = d, byrow = TRUE)
A. <- t(chol(Sigma.))
X.norm. <- rep(mu, each = n) + t(A. %*% t(Z))
plot(rbind(X.norm, X.norm.), xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("black", "royalblue3"), each = n))
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
legend = c(expression(N(mu,Sigma)),
expression("Correlation sign switched")))
## Note: This has all been done based on the same Z (and we also recycle it below)!
### 3 Sampling from the multiv. t distribution (and other normal variance mixtures)
## We recycle the sample with positive correlation
X.norm <- X.norm.
A <- A.
Sigma <- A %*% t(A)
### 3.1 A sample from a multivariate t_nu distribution ##
nu <- 3 # degrees of freedom
set.seed(314)
W <- 1/rgamma(n, shape = nu/2, rate = nu/2) # mixture variable for a multiv. t_nu distribution
plot(W, log = "y")
X.t <- rep(mu, each = n) + sqrt(W) * t(A %*% t(Z))
xlim <- c(-40, 80)
ylim <- c(-20, 60)
plot(rbind(X.t, X.norm), xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("maroon3", "black"), each = n))
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "maroon3"),
legend = c(expression(N(mu,Sigma)),
expression(italic(t)[nu](mu,Sigma))))
## For more about sampling from the multivariate t distribution (especially
## what can go wrong!), see Hofert (2013, "On Sampling from the Multivariate t Distribution")
### 3.2 A 2-point distribution for W (= two 'overlaid' multiv. normals) ##
## Let's change the distribution of W to a 2-point distribution with
## W = 1 with prob. 0.5 and W = 2 with prob. 0.5. What we will see are two 'overlaid'
## normal distributions (as each of W = 1 and W = 2 corresponds to one).
W.binom <- 1 + rbinom(n, size = 1, prob = 0.5)
cols <- rep("royalblue3", n)
cols[W.binom == 2] <- "maroon3"
X.W.binom <- rep(mu, each = n) + sqrt(W.binom) * t(A %*% t(Z))
plot(rbind(X.t, X.W.binom), xlab = expression(X[1]), ylab = expression(X[2]),
col = c(rep("black", n), cols))
legend("topright", bty = "n", pch = rep(1, 3), col = c("royalblue3", "maroon3", "black"),
legend = c(expression(N(mu,W*Sigma)~"for W = 1"),
expression(N(mu,W*Sigma)~"for W = 2"),
expression(italic(t)[nu](mu,Sigma))))
## => With probability 0.5 we sample from a normal variance mixture with W = 1
## and with probability 0.5 we sample from one with W = 2. By using a df for
## W with infinite upper endpoint, we can reach further out in the tails than
## with any multivariate normal distribution by overlaying normals
## with different (unbounded) covariance matrices (and this is what is
## creating heavier tails for the t distribution than the normal).
### 3.3 An *uncorrelated* multivariate t_nu distribution vs an independent one
A. <- diag(sqrt(diag(Sigma))) # Cholesky factor of the diagonal matrix containing the variances
X.t.uncor <- rep(mu, each = n) + sqrt(W) * t(A. %*% t(Z)) # uncorrelated t samples
X.t.ind <- rep(mu, each = n) + matrix(rt(n * d, df = nu), ncol = d) # independent t samples
plot(rbind(X.t.uncor, # uncorrelated t; dependence only introduced by sqrt(W) => shifted spherical distribution
X.t, # full t; dependence introduced by sqrt(W) *and* correlation
X.t.ind), # independence
xlim = xlim, ylim = ylim,
xlab = expression(X[1]), ylab = expression(X[2]),
col = rep(c("royalblue3", "maroon3", "black"), each = n))
legend("topright", bty = "n", pch = rep(1, 3), col = c("black", "royalblue3", "maroon3"),
legend = c(expression("Independent"~italic(t)[nu]),
expression("Uncorrelated"~italic(t)[nu]~"(shifted spherical)"),
expression(italic(t)[nu](mu,Sigma))))
## => The uncorrelated (but not independent) samples show more mass in the joint
## tails than the independent samples. The (neither uncorrelated nor independent)
## t samples show the most mass in the joint tails.
### 3.4 A normal *mean*-variance mixture ##
## Now let's also 'mix' the mean (replace mu by mu(W)), to get a normal
## mean-variance mixture; here: choosing between two different locations mu
## depending on W.
mu. <- t(sapply(W.binom, function(w) mu + if(w == 1) 0 else c(15, 0)))
X.mean.var <- mu. + sqrt(W.binom) * t(A %*% t(Z))
plot(rbind(X.W.binom, X.mean.var), xlab = expression(X[1]), ylab = expression(X[2]),
col = c(rep("black", n), cols))
legend("topleft", bty = "n", pch = rep(1, 3), col = c("black", "royalblue3", "maroon3"),
legend = c(expression(N(mu,W*Sigma)~"for W"%in%"{1,2}"),
"Mean-var. mix for W = 1",
"Mean-var. mix for W = 2"))
## => Not only do the red points show more variation, they *also* have a different
## location; clearly, this sample (blue + red points together) is not
## elliptically distributed anymore.
### 4 Sampling elliptical distributions ##
### 4.1 Sample from a multivariate t_nu via elliptical distributions ##
## Sample from the uniform distribution on the unit sphere (spherical distribution)
set.seed(271)
Z <- matrix(rnorm(n * d), ncol = d)
S <- Z/sqrt(rowSums(Z^2)) # Z/||Z||
lim <- c(-1.3, 1.3)
par(pty = "s")
plot(S, xlim = lim, ylim = lim,
xlab = expression(S[1]), ylab = expression(S[2]))
## Add scale (elliptical distribution) by computing X = RAS (= ARS)
rho <- 0.8 # correlation
P <- matrix(c( 1, rho,
rho, 1), ncol = 2) # correlation matrix
A <- chol(P) # Cholesky factor (here: upper triangular matrix for efficiency reasons)
X. <- t(A %*% t(S)) # (n, d) matrix
par(pty = "s")
plot(X., xlim = lim, ylim = lim,
xlab = expression(X[1]), ylab = expression(X[2]))
## Add radial part of a t_nu distribution (spherical distribution)
nu <- 5
R <- sqrt(d * rf(n, df1 = d, df2 = nu))
X.. <- R * X.
ran <- c(min(X.., -max(X..)), max(X.., -min(X..)))
par(pty = "s")
plot(X.., xlim = ran, ylim = ran, xlab = expression(X[1]), ylab = expression(X[2]))
## Add location (elliptical distribution) by computing X = mu + RAS
mu <- c(ran[2]-max(X..[,1]), ran[2]-max(X..[,2])) # push sample towards the upper right corner
X <- rep(mu, each = n) + X..
par(pty = "s")
plot(X, xlim = ran, ylim = ran, xlab = expression(X[1]), ylab = expression(X[2]))
### 4.2 A 2-point distribution for R ##
R <- 1 + rbinom(n, size = 1, prob = 2/3) # P(R = 1) = 1/3, P(R = 2) = 2/3
cols <- rep("royalblue3", n)
cols[R == 2] <- "maroon3" # prob. 2/3
X <- rep(mu, each = n) + t(A %*% t(R*S)) # same mu, A, S
par(pty = "s")
plot(X, xlab = expression(X[1]), ylab = expression(X[2]), col = cols)
legend("bottomright", bty = "n", pch = c(1,1), col = c("royalblue3", "maroon3"),
legend = c("Elliptical for R = 1", "Elliptical for R = 2"))
## ... in other words, the radial part R scales each point on the ellipse AS
## to lie on another ellipse. For continuously distributed R, these ellipses
## are all overlaid (and thus not visible anymore), but for discrete R, we
## may see them (here: two ellipses).
#### Fitting_a_multivariate_normal_and_t ####
## Fitting a multivariate normal and t distribution to negative log-returns of
## the five components of the Dow Jones index
library(xts) # for time series manipulation
library(QRM) # for fit.mst()
library(mvtnorm) # for sampling from a multivariate normal or t distribution
library(qrmdata) # for Dow Jones constituents data
library(qrmtools) # for returns()
## Load and extract the data we work with and plot
data(DJ_const)
str(DJ_const)
S <- DJ_const['2000-01-01/', c("AAPL", "BA", "INTC", "IBM", "NKE")]
## => We work with Apple, Boeing, Intel, IBM, Nike since 2000-01-01 here
plot.zoo(S, xlab = "Time t")
## Build and plot -log-returns
X <- -returns(S) # compute -log-returns
pairs(as.matrix(X), main = "Scatter plot matrix of risk-factor changes", gap = 0, pch = ".")
plot.zoo(X, xlab = "Time t", main = "")
## For illustration purposes, we treat this data as realizations of iid
## five-dimensional random vectors.
## Fitting a multivariate normal distribution to X and simulating from it
mu <- colMeans(X) # estimated location vector
Sigma <- cov(X) # estimated scale matrix
stopifnot(all.equal(Sigma, var(X)))
P <- cor(X) # estimated correlation matrix
stopifnot(all.equal(P, cov2cor(Sigma)))
n <- nrow(X) # sample size
set.seed(271)
X.norm <- rmvnorm(n, mean = mu, sigma = Sigma) # N(mu, Sigma) samples
## Fitting a multivariate t distribution to X
fit <- fit.mst(X, method = "BFGS") # fit a multivariate t distribution
X.t <- rmvt(n, sigma = as.matrix(fit$Sigma), df = fit$df, delta = fit$mu) # t_nu samples
## Plot (sample from fitted t (red), original sample (black), sample from fitted normal (blue))
dat <- rbind(t = X.t, original = as.matrix(X), norm = X.norm)
cols <- rep(c("maroon3", "black", "royalblue3"), each = n)
pairs(dat, gap = 0, pch = ".", col = cols)
## Pick out one pair (to better see that multivariate t fits better)
plot(dat[,1:2], col = cols) # => the multivariate normal generates too few extreme losses!
legend("bottomright", bty = "n", pch = rep(1, 3), col = c("black", "royalblue3", "maroon3"),
legend = c("-Log-returns", expression("Simulated fitted"~N(mu,Sigma)),
expression("Simulated fitted"~italic(t)[nu](mu,Sigma))))
#### Multivariate_normal ####
## Libraries required
library(QRM)
library(xts)
library(mvtnorm)
library(qrmdata)
library(copula)
## The bivariate normal distribution
BiDensPlot(func = dmnorm, xpt = c(-4, 4), ypts = c(-4, 4),
mu = c(0, 0), Sigma = equicorr(2, -0.7))
## Some trivariate simulated data
rho <- 0.7
P <- matrix(rho,ncol = 3,nrow = 3)+(1-rho)*diag(3)
P
data <- rmvnorm(1000, sigma = P)
## ML estimation
n <- dim(data)[1]
mu.hat <- colMeans(data)
Sigma.hat <- var(data)*(n-1)/n
P.hat <- cor(data)
## Note that var(data) gives unbiased estimator which is not MLE
ll.max <- sum(dmvnorm(data, mean = mu.hat, sigma = Sigma.hat, log = TRUE))
ll.max
## Some 10-d simulated data
rho = 0.6
P <- matrix(rho, ncol = 10, nrow = 10)+(1-rho)*diag(10)
P
data <- rmvnorm(1000, sigma = P)
## Test of marginal normality
apply(data, 2, shapiro.test)
## Tests of joint normality
MardiaTest(data)
jointnormalTest(data)
## Dow Jones data
data("DJ_const")
Sdata <- DJ_const['2000-01-01/',1:10]
Xdata <- diff(log(Sdata))[-1,]
plot.zoo(Xdata)
Xdata.w <- apply.weekly(Xdata,FUN = colSums)
Xdata.m <- apply.monthly(Xdata,FUN = colSums)
Xdata.q <- apply.quarterly(Xdata,FUN = colSums)
## Test of marginal normality
apply(Xdata, 2, shapiro.test)
apply(Xdata.w, 2, shapiro.test)
apply(Xdata.m, 2, shapiro.test)
apply(Xdata.q, 2, shapiro.test)
## Test for joint normality
MardiaTest(as.matrix(Xdata))
jointnormalTest(as.matrix(Xdata))
MardiaTest(as.matrix(Xdata.w))
jointnormalTest(as.matrix(Xdata.w))
MardiaTest(as.matrix(Xdata.m))
jointnormalTest(as.matrix(Xdata.m))
MardiaTest(as.matrix(Xdata.q))
jointnormalTest(as.matrix(Xdata.q))
## Constructing data which are only marginally normal
set.seed(133)
G.cop <- archmCopula("Gumbel", param = 3, dim = 3)
U.G <- rCopula(1000, copula = G.cop)
data = qnorm(U.G)
pairs(data)
apply(data,2,shapiro.test)
MardiaTest(data)
jointnormalTest(data)
#### Multivariate_Student ####
library(xts)
library(mvtnorm)
library(QRM) # for fit.mst()
library(qrmdata)
## Simulate multivariate Student
rho <- 0.7
d <- 3
P <- matrix(rho, ncol = d, nrow = d) + (1-rho)*diag(d)
nu <- 10
data <- rmvt(2000, sigma = P, df = nu)
pairs(data)
## Test of marginal normality
apply(data, 2, shapiro.test)
## Tests of joint normality
MardiaTest(data)
jointnormalTest(data)
## QQplots against Student t
qplotfunc <- function(data){
qqplot(qt(ppoints(data), df = nu), data, xlab = "Theoretical", ylab = "Empirical")
qqline(data, dist = function(p) qt(p, df = nu))
}
par(mfrow = c(1,3))
apply(data, 2, qplotfunc)
par(mfrow = c(1, 1))
## QQplot of Mahalanobis distances (scaled by dimension d)
## against an F(d,nu) distribution
Ddata <- mahalanobis(data, center = FALSE, cov = var(data))/d
qqplot(qf(ppoints(Ddata), df1 = d, df2 = nu), Ddata,
xlab = "Theoretical", ylab = "Empirical")
qqline(Ddata, dist = function(p) qf(p,df1 = d, df2 = nu))
## Estimate multivariate Student
## Dow Jones Data
data("DJ_const")
Sdata <- DJ_const['2000-01-01/',1:10]
Xdata <- diff(log(Sdata))[-1,]
Xdata.w <- apply.weekly(Xdata, FUN = colSums)
Xdata.m <- apply.monthly(Xdata, FUN = colSums)
## Write function for fitting multivariate normal by MLE
fit.mvnorm <- function(data){
n <- dim(data)[1]
mu.hat <- colMeans(data)
Sigma.hat <- var(data)*(n-1)/n
ll.max <- sum(dmvnorm(data, mean = mu.hat, sigma = Sigma.hat, log = TRUE))
list(mu = mu.hat, Sigma = Sigma.hat, ll.max = ll.max)
}
## Fit normal
mod.w.norm = fit.mvnorm(as.matrix(Xdata.w))
mod.w.norm
## Fit Student t
mod.w.t = fit.mst(as.matrix(Xdata.w), method = "BFGS")
mod.w.t
names(mod.w.t)
mod.w.t$df
## Compare
mod.w.t$ll.max - mod.w.norm$ll.max
## Fit normal
mod.m.norm = fit.mvnorm(as.matrix(Xdata.m))
mod.m.norm
## Fit Student t
mod.m.t = fit.mst(as.matrix(Xdata.m), method = "BFGS")
mod.m.t
## Compare
mod.m.t$ll.max - mod.m.norm$ll.max
#### Multivariate_GH ####
## Libraries required
library(xts)
library(qrmdata)
library(ghyp)
## Dow Jones Data
data("DJ_const")
Sdata <- DJ_const['2000-01-01/',1:10]
Xdata <- diff(log(Sdata))[-1,]
Xdata.w <- apply.weekly(Xdata, FUN = colSums) # weekly data
plot.zoo(Xdata.w)
## Fit elliptical models
mod.NIG.sym <- fit.NIGmv(Xdata.w, symmetric = TRUE, nit = 10000)
summary(mod.NIG.sym)
mod.T.sym <- fit.tmv(Xdata.w, symmetric = TRUE, nit = 10000)
summary(mod.T.sym)
mod.GHYP.sym <- fit.ghypmv(Xdata.w, symmetric = TRUE, nit = 10000)
summary(mod.GHYP.sym)
## Fit skewed models
mod.NIG.skew <- fit.NIGmv(Xdata.w, symmetric = FALSE, nit = 10000)
summary(mod.NIG.skew)
mod.T.skew <- fit.tmv(Xdata.w, symmetric = FALSE, nit = 10000)
summary(mod.T.skew)
mod.GHYP.skew <- fit.ghypmv(Xdata.w, symmetric = FALSE, nit = 10000)
summary(mod.GHYP.skew)
mod.NIG.sym@llh
mod.NIG.skew@llh
mod.T.sym@llh
mod.T.skew@llh
mod.GHYP.sym@llh
mod.GHYP.skew@llh
## Hypothesis: elliptical T is good enough
LRstat <- 2*(mod.T.skew@llh-mod.T.sym@llh)
1-pchisq(LRstat,10) # => H0 not rejected
## Hypothesis: elliptical ghyp is good enough
LRstat <- 2*(mod.GHYP.skew@llh-mod.GHYP.sym@llh)
1-pchisq(LRstat,10) # => H0 not rejected
## Hypothesis: skewed t no worse than skewed ghyp
LRstat <- 2*(mod.GHYP.skew@llh-mod.T.skew@llh)
1-pchisq(LRstat,1) # => H0 rejected at 5% level
## Hypothesis: elliptical t no worse than elliptical ghyp
LRstat <- 2*(mod.GHYP.sym@llh-mod.T.sym@llh)
1-pchisq(LRstat,1) # => H0 rejected at 5% level
## Find skewness parameters
getSlots("mle.ghyp")
mod.GHYP.skew@gamma
## Make picture of two returns
Xdata <- Xdata.w[,c(5,8)]
plot(as.numeric(Xdata[,1]), as.numeric(Xdata[,2]),
xlab = names(Xdata)[1], ylab = names(Xdata)[2])
mod <- fit.ghypmv(Xdata, symmetric = FALSE, nit = 10000)
summary(mod)
xvals <- seq(from = min(Xdata[,1]), to = max(Xdata[,2]), length.out = 50)
yvals <- seq(from = min(Xdata[,2]), to = max(Xdata[,2]), length.out = 50)
outer.func <- function(x,y) dghyp(cbind(x,y), mod)
zvals <- outer(xvals, yvals, outer.func)
zvals.un <- unique(as.numeric(zvals))
contour(xvals, yvals, zvals, add = TRUE, col = 2,
levels = seq(from = quantile(zvals.un,0.6), to = quantile(zvals.un,0.9), length = 10))
#### Fundamental_factor_model ####
## Fitting a fundamental factor model using industry sector information
library(xts) # for time series manipulation
library(qrmdata) # for SP500 (constituents) data
library(qrmtools) # for plot_NA() and plot_matrix()
### 1 Data preparation ##
## Load and extract the data we work with (all available since 1995)
## Index
data(SP500_const) # index data
SP500.const <- SP500_const['1995-01-01/',] # all since 1995
dim(SP500.const)
NA_plot(SP500.const[,1:100]) # => many stocks have missing values (look at first 100)
SP500.const <- SP500.const[, colSums(is.na(SP500.const)) <= 0.1 * nrow(SP500.const)] # omit columns with more than 10% NA
SP500.const <- na.fill(SP500.const, fill = "extend") # fill the remaining NAs
dim(SP500.const)
## Caution: only execute next block if you want to see all the data plotted!!!
par(ask = TRUE)
nplots <- ceiling(dim(SP500.const)[2]/20)
for (i in 1:nplots){
r1 <- 20*(i-1)+1
r2 <- min(20*i,dim(SP500.const)[2])
plot.zoo(SP500.const[,r1:r2], xlab = "Time t", main = "SP500 Constituents")
}
par(ask = FALSE)
## Remove certain stocks
## AIG, C, ETFC (badly affected by financial crisis)
## T, CTL, FTR, VZ (only 4 telecommunications companies so this sector is thin)
## CMCSK (moves together with CMCSA)
## FOX (moves together with FOXA)
ignore <- list("AIG","T","CTL","C","CMCSK","ETFC","FTR","FOX","VZ")
which(names(SP500.const) %in% ignore)
SP500.const <- SP500.const[,-which(names(SP500.const) %in% ignore)]
dim(SP500.const)
## Build and plot log-returns
X.const <- diff(log(SP500.const))[-1,] # compute -log-returns
## We use monthly data here as basis (and compute monthly log-returns)
X <- apply.monthly(X.const, FUN = colSums) # 252 by 367 matrix
plot.zoo(X[,1:20], type = "h", xlab = "Time t",
main = "Monthly risk-factor changes (log-returns) of SP500 constituents")
summary(SP500_const_info)
firms <- as.character(SP500_const_info$Ticker)
## Make two firm names consistent
firms[firms == "BRK-B"] <- "BRK.B"
firms[firms == "BF-B"] <- "BF.B"
sectors <- as.character(SP500_const_info$Sector)
## Construct a factor X.sectors and tabulate levels
X.sectors <- sectors[which(firms %in% colnames(X))]
X.sectors <- as.factor(X.sectors)
table(X.sectors)
levels(X.sectors) <- c("Con-Disc.","Con-Stap.","Energy",
"Financials","Health","Industrials","IT",
"Materials","Utilities")
### 2 Model fitting ##
## Fit a cross-sectional regression model X_t = B*F_t + eps
## A regression is fitted at each time point
## No intercept is required (-1)
mod <- lm(t(X) ~ X.sectors -1)
B <- model.matrix(mod)
## Factors can be constructed from coefficient matrix
coef.mat <- t(coef(mod))
dimnames(coef.mat)[[2]] <- levels(X.sectors)
F <- xts(coef.mat,time(X))
plot.zoo(F,type = "h")
## Calculate error variances for each stock
eps <- t(residuals(mod))
eps.variances <- diag(var(eps))
## Display them from smallest to largest
sort(eps.variances)
## Use generalized least squares to obtain better estimates
## See documentation of lm for more details
mod2 <- lm(t(X) ~ X.sectors -1, weights = 1/eps.variances)
coef.mat <- t(coef(mod2))
dimnames(coef.mat)[[2]] <- levels(X.sectors)
F <- xts(coef.mat,time(X))
plot.zoo(F,type = "h")
cor(F)
eps <- t(residuals(mod))
## This matrix should be near-diagonal for perfect factor model
## Of course this won't be the case in practice
cor.eps <- cor(eps)
## Is Cor(eps) (roughly) diagonal? Check upper corner
diag(cor.eps) <- NA
matrix_plot(cor.eps[1:75,1:75])
## Are the errors uncorrelated with the factors?
cor.eps.F <- cor(eps, F)
summary(cor.eps.F) # => there are some large correlations
## Construct the implied covariance and correlation matrix
Ups <- cov(eps) # Upsilon (covariance matrix of epsilon)
Omega <- as.matrix(cov(F)) # Omega (covariance matrix of F; systematic risk)
Sigma <- B %*% Omega %*% t(B) + diag(diag(Ups)) # Cov(X)
P <- cov2cor(Sigma) # Cor(X)
## Look at discrepancies between the factor model correlation matrix and the
## sample correlation matrix; check upper corner
err <- P-cor(X)
matrix_plot(err[1:75,1:75])
#### Macroeconomic_factor_model ####
## Fitting a factor model
library(xts) # for time series manipulation
library(qrmdata) # for Dow Jones (constituents) data
library(qrmtools) # for plot_NA() and plot_matrix()
### 1 Data preparation ##
## Load and extract the data we work with (all available since 1990) and plot
## Index
data(DJ) # index data
DJ. <- DJ['1990-01-01/'] # all since 1990
plot.zoo(DJ., xlab = "Time t", ylab = "Dow Jones Index")
## Constituents
data(DJ_const) # constituents data
DJ.const <- DJ_const['1990-01-01/',] # all since 1990
NA_plot(DJ.const) # => use all but the two columns with lots of NAs
DJ.const <- DJ.const[, colSums(is.na(DJ.const)) <= 0.1 * nrow(DJ.const)] # omit columns with more than 10% NA
DJ.const <- na.fill(DJ.const, fill = "extend") # fill the remaining NAs
plot.zoo(DJ.const, xlab = "Time t", main = "Dow Jones Constituents")
## Build and plot log-returns
## Index
X. <- diff(log(DJ.))[-1,] # compute -log-returns
plot.zoo(X., xlab = "Time t", ylab = expression(X[t]), main = "Risk-factor changes (log-returns) of Dow Jones index")
## Constituents
X.const <- diff(log(DJ.const))[-1,] # compute -log-returns
if(FALSE) # more time-consuming
pairs(as.matrix(X.const), gap = 0, pch = ".",
main = "Scatter plot matrix of risk-factor changes (log-returns) of Dow Jones constituents")
plot.zoo(X.const, xlab = "Time t", main = "Risk-factor changes (log-returns) of Dow Jones constituents")
## We use monthly data here as basis (and compute monthly log-returns)
X <- apply.monthly(X.const, FUN = colSums) # (312, 28)-matrix
plot.zoo(X, type = "h", xlab = "Time t", main = "Monthly risk-factor changes (log-returns) of Dow Jones constituents")
F <- apply.monthly(X., FUN = colSums)
plot(F, type = "h", xlab = "Time t", ylab = expression(X[t]), main = "Monthly risk-factor changes (log-returns) of Dow Jones index")
### 2 Model fitting ##
## Fit a multivariate regression model X = a + B*F + eps
(res <- lm(X ~ F)) # more details via summary()
## Get parameter estimates
names(res)
par.ests <- coefficients(res)
a <- par.ests[1,]
B <- par.ests[2,]
# tail(X[,1:3])
# tail(predict(res)[,1:3])
# forecast::forecast(res, newdata = F[length(F)], h = 1)
## Get residuals and estimate their correlation matrix
eps <- resid(res) # (312, 28)-matrix
cor.eps <- cor(eps)
## Is Cor(eps) (roughly) diagonal?
matrix_plot(cor.eps, at = seq(-1, 1, length.out = 200)) # => yes (as required)
## Are the errors uncorrelated with the factors?
cor.eps.F <- cor(eps, F) # 28 correlations (idiosyncratic risk)
summary(cor.eps.F) # => yes (as required)
## Construct the implied covariance and correlation matrix
Ups <- cov(eps) # Upsilon (covariance matrix of epsilon)
Omega <- as.matrix(cov(F)) # Omega (covariance matrix of F; systematic risk)
Sigma <- B %*% Omega %*% t(B) + diag(diag(Ups)) # Cov(X); (28, 28)-matrix
rownames(Sigma) <- colnames(Sigma)
P <- cov2cor(Sigma) # Cor(X)
## Look at discrepancies between the factor model correlation matrix and the
## sample correlation matrix
err <- P-cor(X)
matrix_plot(err, at = seq(-1, 1, length.out = 200))
#### PCA_factor_model ####
## Fitting a PCA factor model
library(xts) # for time series manipulation
library(qrmdata) # for Dow Jones (constituents) data
library(qrmtools) # for plot_NA() and plot_matrix()
### 1 Data preparation ###
## Load and extract the data we work with (all available since 1990) and plot
## Index
data(DJ) # index data
DJ. <- DJ['1990-01-01/'] # all since 1990
plot.zoo(DJ., xlab = "Time t", ylab = "Dow Jones Index")
## Constituents
data(DJ_const) # constituents data
DJ.const <- DJ_const['1990-01-01/',] # all since 1990
NA_plot(DJ.const) # => use all but the two columns with lots of NAs
DJ.const <- DJ.const[, colSums(is.na(DJ.const)) <= 0.1 * nrow(DJ.const)] # omit columns with more than 10% NA
DJ.const <- na.fill(DJ.const, fill = "extend") # fill the remaining NAs
plot.zoo(DJ.const, xlab = "Time t", main = "Dow Jones Constituents")
## Build and plot log-returns
## Index
X. <- diff(log(DJ.))[-1,] # compute -log-returns
plot.zoo(X., xlab = "Time t", ylab = expression(X[t]),
main = "Risk-factor changes (log-returns) of Dow Jones index")
## Constituents
X.const <- diff(log(DJ.const))[-1,] # compute -log-returns
if(FALSE) # more time-consuming
pairs(as.matrix(X.const), gap = 0, pch = ".",
main = "Scatter plot matrix of risk-factor changes (log-returns) of Dow Jones constituents")
plot.zoo(X.const, xlab = "Time t", main = "Risk-factor changes (log-returns) of Dow Jones constituents")
## We use monthly data here as basis (and compute monthly log-returns)
X <- apply.monthly(X.const, FUN = colSums) # (312, 28)-matrix
plot.zoo(X, type = "h", xlab = "Time t",
main = "Monthly risk-factor changes (log-returns) of Dow Jones constituents")
Xindex <- apply.monthly(X., FUN = colSums)
plot(Xindex, type = "h", xlab = "Time t", ylab = expression(X[t]),
main = "Monthly risk-factor changes (log-returns) of Dow Jones index")
### 2 Model fitting ###
## Principal component analysis
PCA <- princomp(X) # determine the principal components
PCA. <- prcomp(X) # another option
summary(PCA)
plot(PCA) # show important components
loadings(PCA) # factor loadings (representing how much a factor explains a variable)
## Working with the principal components
mu <- PCA$center # estimated centers
Y <- PCA$scores # estimated principal components of X
nprin <- 3 # number of important principal components
Y1 <- xts(Y[,1:nprin], time(X)) # grab out the first so-many important principal components of X
plot.zoo(Y1, type = "h", xlab = "Time", ylab = paste("Component", 1:nprin),
main = "Principal components of X")
cor(Y1)
Y2 <- Y[,(nprin+1):ncol(Y)] # ignored principal components of X
## Interpreting first principal component
plot.zoo(merge(Xindex,Y1[,1]),main=" ") # compare DJ index returns and first PC
plot(as.numeric(Xindex),-as.numeric(Y1[,1]),xlab="Index",ylab="PC1 (neg)")
## Reconstruct X from Y1 and Y2
G <- unclass(loadings(PCA)) # compute G, see McNeil, Frey, Embrechts (2015, (6.64))
G1 <- G[,1:nprin] # compute G_1, see McNeil, Frey, Embrechts (2015, (6.65))
G2 <- G[,(nprin+1):ncol(G)]
eps <- G2 %*% t(Y2) # compute epsilon
eps.cor <- cor(t(eps)) # compute correlations of epsilons
diag(eps.cor) <- NA # neglect diagonal elements
matrix_plot(eps.cor) # not perfectly diagonal
X. <- t(mu + G1 %*% t(Y1) + eps)
err <- X.-X
head(err) # check that differences are zero vectors
summary(as.vector(err))
## Note: We could fit a univariate GARCH model to each PCA factor.
## This is called PC-GARCH model.
#### Projections_of_uniform_on_sphere ####
## Orthogonal projections of the uniform distribution on the unit sphere S_{d-1}
## in d dimensions down by two dimensions are uniform in the (d-2)-ball;
## see http://www.math.univ-toulouse.fr/~letac/Archimede.pdf
## Let's check if that's the case!
library(rgl)
## Sample from the uniform distribution on the unit d-sphere S_{d-1}
n <- 5e3 # sample size
d <- 4 # maximal dimension investigated
set.seed(271) # for reproducibility
Z <- matrix(rnorm(n*d), ncol = d) # generate independent N(0,1)
### 1 Projections from d = 3 to d = 2 ##
S2 <- Z[,1:3]/sqrt(rowSums(Z[,1:3]^2)) # ~ U(S_2)
pairs(S2, gap = 0, pch = ".") # => does not look uniform (fine)
### 2 Projections from d = 3 to d = 1 ##
plot(S2[,1], ylab = expression(S[1]), cex = 0.5) # => looks uniform in [0,1]
plot(S2[,2], ylab = expression(S[2]), cex = 0.5) # => looks uniform in [0,1]
plot(S2[,3], ylab = expression(S[3]), cex = 0.5) # => looks uniform in [0,1]
### 3 Projections from d = 4 to d = 3 ##
S3 <- Z[,1:4]/sqrt(rowSums(Z[,1:4]^2)) # ~ U(S_3)
colnames(S3) <- paste0("S", 1:4)
plot3d(S3[,c(2,3,4)]) # => looks uniform in the ball in 3d but difficult to judge
S3. <- S3[0.2 < S3[,2] & S3[,2] <= 0.3, c(2,3,4)] # pick out a slice (0.2, 0.3]
plot(S3.) # slice looks uniform
plot3d(S3[,c(1,3,4)])
plot3d(S3[,c(1,2,4)])
plot3d(S3[,c(1,2,3)])
## Note: Although not covered by the above result, it *seems* that projections
## from U(S_3) to IR^3 are also uniform in the 3-ball. More investigations
## necessary.
### 4 Projections from d = 4 to d = 2 ##########################################
pairs(S3, gap = 0, pch = ".", # => looks (pairwise) uniform in the unit disk
labels = as.expression(sapply(1:4, function(j) bquote(S[.(j)]))))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.