MAT 4376F/5313F
R – Fitting multivariate distributions
library(mvtnorm);
library(MASS);
library(QRM);
library(evir);
# Simulate bivariate normal data
mu <- c(0,0) # Mean
Sigma <- matrix(c(1, .9, .9, 1), 2) # Covariance matrix
n=1000;
Xnorm<- rmvnorm(n, mean=mu, sigma=Sigma, method="chol")
par(mfrow=c(1,2))
qqnorm(Xnorm[,1]); qqnorm(Xnorm[,2]); # Marginally data are normal
par(mfrow=c(1,2))
hist(Xnorm[,1]); hist(Xnorm[,2]);
# scatter plot
par(mfrow=c(1,1))
plot(Xnorm[,1],Xnorm[,2]);
# Calculate kernel density estimate
Bivnorm.kde <- kde2d(Xnorm[,1], Xnorm[,2], n = 50) # from MASS package
par(mfrow=c(1,1))
persp(Bivnorm.kde, phi = 45, theta = 30, shade = .1, border = NA);
# Contour plot for bivariate normal density
x.points <-seq(-3,3,length.out=100);
y.points <- x.points;
z <-matrix(0,nrow=100,ncol=100);
mu <-c(0,0);
for (i in 1:100)
{for (j in 1:100)
{z[i,j] <-dmvnorm(c(x.points[i],y.points[j]),mean=mu,sigma=Sigma)}}
contour(x.points,y.points,z)
# Contour plot for the bivariate data
contour(Bivnorm.kde)
# Simulate bivariate t
Sigma <- matrix(c(1, .9, .9, 1), 2) # Sigma matrix
nu <- 4 # number of degrees of freedom
Xt<- rmvt(n, sigma=Sigma, df=nu);
par(mfrow=c(1,2))
qqnorm(Xt[,1]); qqnorm(Xt[,2]); # Marginally data are not normal
# scatter plot
par(mfrow=c(1,1))
plot(Xt[,1],Xt[,2]);
# Calculate kernel density estimate
Bivt.kde <- kde2d(Xt[,1], Xt[,2], n = 50) # from MASS package
par(mfrow=c(1,1))
persp(Bivt.kde, phi = 45, theta = 30, shade = .1, border = NA);
# Comparison of bivariate normal and t
par(mfrow=c(1,2))
plot(Xnorm[,1],Xnorm[,2]);
plot(Xt[,1],Xt[,2]);
par(mfrow=c(1,2))
persp(Bivnorm.kde, phi = 45, theta = 30, shade = .1, border = NA);
persp(Bivt.kde, phi = 45, theta = 30, shade = .1, border = NA);
par(mfrow=c(1,2))
contour(Bivnorm.kde)
contour(Bivt.kde)
# Fitting bivariate t distribution
fit <- fit.mst(Xt) # fit a multivariate t distribution
mu.fit <- fit$mu # estimated location vector;
Sigma,fit <- as.matrix(fit$Sigma) # estimated scale matrix;
nu.fit <- fit$df # estimated degrees of freedom
fit <- fit.mst(Xnorm) # fit a multivariate t distribution to normal data
mu.fit <- fit$mu # estimated location vector;
Sigma.fit <- as.matrix(fit$Sigma) # estimated scale matrix;
nu.fit <- fit$df # estimated degrees of freedom
# Working with real data
data(bmw); data(siemens);
X1=-bmw; X2=-siemens;
par(mfrow=c(1,2));
plot.ts(X1); plot.ts(X2);
par(mfrow=c(1,2));
qqnorm(X1); qqnorm(X2);
par(mfrow=c(1,1));
plot(X1,X2);
Biv.kde <- kde2d(X1, X2, n = 50) # from MASS package
par(mfrow=c(1,1))
persp(Biv.kde, phi = 45, theta = 30, shade = .1, border = NA);
par(mfrow=c(1,1))
contour(Biv.kde)
X=cbind(X1,X2);
fit <- fit.mst(X) # fit a multivariate t distribution
mu.fit <- fit$mu # estimated location vector;
Sigma.fit <- as.matrix(fit$Sigma) # estimated scale matrix;
nu.fit <- fit$df # estimated degrees of freedom
# "Goodness of fit test"
n=length(X1);
Xnorm<- rmvnorm(n, mean=mu.fit, sigma=Sigma.fit, method="chol")
Xt<- rmvt(n, sigma=Sigma.fit, df=nu.fit);
Bivnorm.kde <- kde2d(Xnorm[,1], Xnorm[,2], n = 50);
Bivt.kde <- kde2d(Xt[,1], Xt[,2], n = 50);
Biv.kde <- kde2d(X1, X2, n = 50);
par(mfrow=c(1,3));
persp(Biv.kde, phi = 45, theta = 30, shade = .1, border = NA);
persp(Bivnorm.kde, phi = 45, theta = 30, shade = .1, border = NA);
persp(Bivt.kde, phi = 45, theta = 30, shade = .1, border = NA);
par(mfrow=c(1,3));
contour(Biv.kde)
contour(Bivnorm.kde)
contour(Bivt.kde)