# Joint Distribution Of Stock Returns [Part I]

While linear correlation represents the appropriate measure of dependence for normally distributed asset returns, capturing co-movements between non-normal financial asset returns requires the theory of copulas, a flexible way for modelling general multivariate dependencies. Let’s download 5 years worth of price data for JPM and MS, set up the data environment, calculate log returns and explore the marginal distribution of returns.

```library(fExtremes)
library(timeSeries)
library(LambertW)
library(quantmod)
library(MSBVAR)
library(spd)
library(PerformanceAnalytics)

ticker<-c('JPM','MS')
symbols<-unlist(strsplit(ticker,split=","))
getSymbols(symbols,auto.assign=TRUE,from='1996-01-01',to='2011-12-31')
ret.mat<-as.timeSeries(cbind(tail(diff(log(get(symbols[1])[,4])),-1),tail(diff(log(get(symbols[2])[,4])),-1)))
ret<-paste(ticker,"Returns")
ret.unl<-unlist(strsplit(ret,","))
colnames(ret.mat)<-ret.unl

# Ensure return series contains only non zero values
nonz.ind<-(seriesData(ret.mat[,1]) !=0 & seriesData(ret.mat[,2])!=0)
ret.mat.nz<-ret.mat[nonz.ind,]

# ##############################################################################
# Exploring the marginal distribution of returns
# ##############################################################################

# [1]Plotting the time series

windows()
layout(c(1,1,2,2))
par(mai=c(0.3,0.75,0,0.1))
plot(ret.mat.nz[,1],cex=0.7,cex.axis=0.9,cex.lab=1,cex.main=0.7,col='steelblue',font.axis=2,xaxt='n',xlab='',bty='n')
abline(h=0)
par(mai=c(0.5,0.75,0.3,0.1))
plot(ret.mat.nz[,2],cex=0.7,cex.axis=0.9,cex.lab=1,cex.main=0.7,col='darkblue',font.axis=2,xaxt='n',xlab='',bty='n')
abline(h=0)

# [2] Normality tests
windows()
x_min<-min(ret.mat.nz)
x_max<-max(ret.mat.nz)
layout(c(1,1,2,2,3,3))
par(mai=c(0.1,0.75,0.3,0.1))
par(mai=c(0.1,0.75,0,0.1))
par(mai=c(0.3,0,0,0.1))

windows()
chart.QQPlot(ret.mat.nz[,1],distribution='norm',line="quartiles",element.color="black",cex.lab=0.75,cex.axis=0.75,cex.main=0.75,lwd=1,pch=1,cex=0.5)

windows()
chart.QQPlot(ret.mat.nz[,2],distribution='norm',line="quartiles",element.color="black",cex.lab=0.75,cex.axis=0.75,cex.main=0.75,lwd=1,pch=1,cex=0.5)```

Created by Pretty R at inside-R.org

Visual inspection of the above figures should convince us that neither JPM nor MS returns are normally distributed, invalidating the use of linear correlation as a reliable measure of dependency between these financial assets. In addition to the univariate marginal distribution explored above, let’s look at the joint bivariate distribution of MS and JPM log-returns and compare it to simulated bi-variate normal data.

```# [3] Compare Empirical Scatter with simulated bivariate normal data

simulated.ret<-rmultnorm(n=nrow(ret.mat.nz),vmat=var(ret.mat.nz),mu=colMeans(ret.mat.nz))
colnames(simulated.ret)<-paste("Simulated Bivariate Normal Return for",colnames(ret.mat.nz))
tot.ret<-cbind(ret.mat.nz,simulated.ret)
emp.corr<-cor(ret.mat.nz[,1],ret.mat.nz[,2])

windows(width=7*2)
layout(matrix(c(1,2,3,4,5,6),nrow=2,ncol=3))
par(mai=c(0.7,0.7,0.3,0.1))
plot(ylim=c(x_min,x_max),xlim=c(x_min,x_max),xlab=ret.unl[1],ylab=ret.unl[2],x=ret.mat.nz[,1],y=ret.mat.nz[,2],pch=1,col="darkblue",cex.main=0.85,cex.lab=0.85,cex=0.85,cex.axis=0.85,main=paste("Empirical data with rho of\n ",emp.corr))
grid()
abline(h=0,v=0)

par(mai=c(0.7,0.7,0.3,0.1))
plot(ylim=c(x_min,x_max),xlim=c(x_min,x_max),xlab=paste("Simulated ",ret.unl[1]),ylab=paste("Simulated ",ret.unl[2]),x=simulated.ret[,1],y=simulated.ret[,2],pch=1,col="darkblue",cex.main=0.85,cex.lab=0.75,cex=0.85,cex.axis=0.85,main="Simulated data from \n Bivariate Normal Distribution")
grid()
abline(h=0,v=0)

par(mai=c(0.7,0.7,0.3,0.1))

par(mai=c(0.7,0.7,0.3,0.1))

par(mai=c(0.7,0.7,0.3,0.1))

par(mai=c(0.7,0.7,0.3,0.1))

Created by Pretty R at inside-R.org

Visual inspection of the graphs suggest that the simulated bi-variate normal data matches the joint behaviour of returns in the middle of the distribution quite well but does not capture the observed dependence in the tails. When the bi-variate normal distribution does not adequately describe joint returns,pearson correlation may not be a proper measure of dependence. Tails of asset returns can be modelled using the generalised pareto distribution(parametric) while the centre of distribution can be captured using an Empirical distribution function (non parametric). This results in a semi-parametric approach to modelling the bi-variate data. Incidentally, if we had imposed a normality assumption on the JPM and MS data, data points would be clustered around 0 in the expected way (last 3 graphs).

Let’s fit a semi parametric model to the data and compare the actual returns with simulated normal, and semi parametric  returns. The semi parametric model is a better fit than the normal model.

```# [4] Modelling the Tails of the distribution

# [4.a] The Mean Excess Plot should tell us where to cut the tails from the entire distribution
windows()
layout(matrix(c(1,2),ncol=1,nrow=2))
par(mai=c(1,0.75,0.3,0.3))
mePlot(ret.mat.nz[,1],cex=0.75,cex.axis=0.75,cex.main=0.85,cex.lab=0.75,sub=colnames(ret.mat.nz)[1])
mePlot(ret.mat.nz[,2],cex=0.75,cex.axis=0.75,cex.main=0.85,cex.lab=0.75,sub=colnames(ret.mat.nz)[2])

# [4.b] Perform a sermiparametric fit; Normal density for the middle of the distribution, GPD for the tails

semiparam.1.fit<-spdfit(ret.mat.nz[,1],upper=0.9,lower=0.1,tailfit='GPD',kernelfit='normal')
semiparam.2.fit<-spdfit(ret.mat.nz[,2],upper=0.9,lower=0.1,tailfit='GPD',kernelfit='normal')

# Plots for the first series
windows()
layout(matrix(c(1,2),nrow=2,ncol=1))
par(mai=c(0.75,0.75,0.7,0.3),cex.main=0.85,cex.lab=0.85,cex=0.75,cex.axis=0.85,bty='n')
plot(semiparam.1.fit,which=1)
mtext(paste(ticker[1],"Data"),side=3,col="darkblue",cex=0.75)
plot(semiparam.1.fit,which=2)
mtext(paste(ticker[1],"Data"),side=3,col="darkblue",cex=0.75)

windows()
layout(matrix(c(1,2),nrow=2,ncol=1))
par(mai=c(0.5,0.75,0.7,0.3),cex.main=0.85,cex.lab=0.85,cex=0.75,cex.axis=0.85,bty='n')
plot(semiparam.1.fit,which=3)
mtext(paste("Fitting the Upper Tails of the ",ticker[1],"Return Series"),side=1,col="darkblue",cex=1,outer=T,font=4)

windows()
layout(matrix(c(1,2),nrow=2,ncol=1))
par(mai=c(0.75,0.75,0.7,0.3),cex.main=1,cex.lab=0.85,cex=0.75,cex.axis=0.85,bty='n')
plot(semiparam.1.fit,which=4)
mtext(paste("Fitting the Lower Tails of the ",ticker[1],"Return Series"),side=1,col="darkblue",cex=1,outer=T,font=4)

# Plots for the second Series

windows()
layout(matrix(c(1,2),nrow=2,ncol=1))
par(mai=c(0.75,0.75,0.7,0.3),cex.main=0.85,cex.lab=0.85,cex=0.75,cex.axis=0.85,bty='n')
plot(semiparam.2.fit,which=1)
mtext(paste(ticker[2],"Data"),side=3,col="darkred",cex=0.75)
plot(semiparam.2.fit,which=2)
mtext(paste(ticker[2],"Data"),side=3,col="darkred",cex=0.75)

windows()
layout(matrix(c(1,2),nrow=2,ncol=1))
par(mai=c(0.5,0.75,0.7,0.3),cex.main=0.85,cex.lab=0.85,cex=0.75,cex.axis=0.85,bty='n')
plot(semiparam.2.fit,which=3)
mtext(paste("Fitting the Lower Tails of the ",ticker[2],"Return Series"),side=1,col="darkred",cex=1,outer=T,font=4)

windows()
layout(matrix(c(1,2),nrow=2,ncol=1))
par(mai=c(0.5,0.75,0.7,0.3),cex.main=0.85,cex.lab=0.85,cex=0.75,cex.axis=0.85,bty='n')
plot(semiparam.2.fit,which=4)
mtext(paste("Fitting the Lower Tails of the ",ticker[2],"Return Series"),side=1,col="darkred",cex=1,outer=T,font=4)

#Comparing the scatterplots for
# actual returns,simulated bivariate normal returns,and simulated Semiparametric returns

semiparam.1.sim<-rspd(nrow(ret.mat.nz),semiparam.1.fit)
semiparam.2.sim<-rspd(nrow(ret.mat.nz),semiparam.2.fit)
semiparametric.sim.ret<-cbind(semiparam.1.sim,semiparam.2.sim)
colnames(semiparametric.sim.ret)<-paste("Semiparametric simulated returns for",ticker)

windows(width=7*2,height=5)
layout(matrix(c(1,2,3),nrow=1,ncol=3))
plot(ylim=c(x_min,x_max),xlim=c(x_min,x_max),xlab=colnames(ret.mat.nz)[1],ylab=colnames(ret.mat.nz)[2],cex.main=1,cex=1,cex.lab=0.95,cex.axis=0.85,x=ret.mat.nz[,1],y=ret.mat.nz[,2],main=paste("Actual Return Data For \n",ticker[1],"and",ticker[2]),col='blue',pch=10)
abline(h=0,v=0)
grid()

plot(ylim=c(x_min,x_max),xlim=c(x_min,x_max),xlab=colnames(simulated.ret)[1],ylab=colnames(simulated.ret)[2],cex.main=1,cex=1,cex.lab=0.95,cex.axis=0.85,x=simulated.ret[,1],y=simulated.ret[,2],main=paste("Simulated Return Data From \n","BiVariate Normal Distribution"),col='darkblue',pch=10)
abline(h=0,v=0)
grid()

plot(ylim=c(x_min,x_max),xlim=c(x_min,x_max),xlab=colnames(semiparametric.sim.ret)[1],ylab=colnames(semiparametric.sim.ret)[2],cex.main=1,cex=1,cex.lab=0.95,cex.axis=0.85,x=semiparametric.sim.ret[,1],y=semiparametric.sim.ret[,2],main=paste("Simulated Return Data From \n","Semiparametric GPD"),col='darkblue',pch=10)
abline(h=0,v=0)
grid()```

Created by Pretty R at inside-R.org