---
title: "Descriptive Statistics for Financial Time Series"
author: "Eric Zivot"
date: "Monday, April 27, 2015"
output: slidy_presentation
---

## Set options and load packages

```{r}
options(digits=4, width=70)
library(car)
library(corrplot)
library(dygraphs)
library(IntroCompFinR)
library(PerformanceAnalytics)
library(sn)
library(tseries)
library(zoo)
Sys.setenv(TZ="UTC")
```

## Get data from introCompfinR package

```{r}
# get daily prices on Microsoft and S&P500 Index from 1993-01-03 to 2014-12-31
data(msftDailyPrices, sp500DailyPrices)
str(msftDailyPrices)
# set sample for descriptive statistical analysis
smpl = "1998-01::2012-05"
msftDailyPrices = msftDailyPrices[smpl]
sp500DailyPrices = sp500DailyPrices[smpl]
# extract end-of-month prices
msftMonthlyPrices = to.monthly(msftDailyPrices, OHLC=FALSE)
sp500MonthlyPrices = to.monthly(sp500DailyPrices, OHLC=FALSE)
```

## Create xts objects with both series

```{r}
msftSp500DailyPrices = merge(msftDailyPrices, sp500DailyPrices)
msftSp500MonthlyPrices = merge(msftMonthlyPrices, sp500MonthlyPrices)
# show data
head(msftSp500DailyPrices, n=3)
head(msftSp500MonthlyPrices, n=3)

```

## Calculate simple and cc returns

```{r}
msftMonthlyRetS = Return.calculate(msftMonthlyPrices, method="simple")
msftDailyRetS = Return.calculate(msftDailyPrices, method="simple")
sp500MonthlyRetS = Return.calculate(sp500MonthlyPrices, method="simple")
sp500DailyRetS = Return.calculate(sp500DailyPrices, method="simple")
msftSp500MonthlyRetS = Return.calculate(msftSp500MonthlyPrices, method="simple")
msftSp500DailyRetS = Return.calculate(msftSp500DailyPrices, method="simple")
# show data
head(msftSp500DailyRetS, n=3)
head(msftSp500MonthlyRetS, n=3)

```

## Remove first NA observations

```{r}
# Note: can alternatively use na.omit() to eliminate NA values
msftMonthlyRetS = msftMonthlyRetS[-1]
msftDailyRetS = msftDailyRetS[-1]
sp500MonthlyRetS = sp500MonthlyRetS[-1]
sp500DailyRetS = sp500DailyRetS[-1]
msftSp500MonthlyRetS = msftSp500MonthlyRetS[-1]
msftSp500DailyRetS = msftSp500DailyRetS[-1]

```

## Compute cc returns

```{r}
msftMonthlyRetC = log(1 + msftMonthlyRetS)
sp500MonthlyRetC = log(1 + sp500MonthlyRetS)
msftSp500MonthlyRetC = merge(msftMonthlyRetC, sp500MonthlyRetC)
# show data
head(msftSp500MonthlyRetC, n=3)

```

## Plot prices in two panel graph

```{r}
plot.zoo(msftSp500MonthlyPrices, main="", lwd=2, col="blue")

```

## Plot returns in two panel graph

```{r}
my.panel <- function(...) {
  lines(...)
  abline(h=0)
}
plot.zoo(msftSp500MonthlyRetS, main="", panel=my.panel, 
         lwd=2, col="blue")

```

## Plot returns in one panel and add a horizontal line

```{r}
plot.zoo(msftSp500MonthlyRetS, plot.type="single", main="",
     col = c("red", "blue"), lty=c("dashed", "solid"), 
     lwd=2, ylab="Returns")
abline(h=0)
legend(x="bottomright", legend=colnames(msftSp500MonthlyRetS), 
       lty=c("dashed", "solid"), lwd=2, 
       col=c("red","blue"))

```

## Comparing simple and cc returns

```{r}
retDiff = msftMonthlyRetS - msftMonthlyRetC
dataToPlot = merge(msftMonthlyRetS, msftMonthlyRetC, retDiff)
plot.zoo(dataToPlot, plot.type="multiple", main="",
     panel = my.panel,
     lwd=2, col=c("black", "blue", "red"))
```

## Plotting daily returns

```{r}
plot.zoo(msftSp500DailyRetS, main="",
     panel=my.panel, col=c("black", "blue"))
```

## Equity curves for MSFT and SP500

```{r}
equityCurveMsft = cumprod(1 + msftMonthlyRetS)  
equityCurveSp500 = cumprod(1 + sp500MonthlyRetS)  
dataToPlot = merge(equityCurveMsft, equityCurveSp500)
plot.zoo(dataToPlot, plot.type="single", ylab="Cumulative Returns",
     col=c("black", "blue"), lwd=2)
legend(x="topleft", legend=c("MSFT", "SP500"),
       col=c("black", "blue"), lwd=2)
```

## Histograms of daily and monthly returns

```{r}
par(mfrow=c(2,2))
  hist(msftMonthlyRetS, main="", col="cornflowerblue")
  hist(msftDailyRetS, main="", col="cornflowerblue")
  hist(sp500MonthlyRetS, main="", col="cornflowerblue")
  hist(sp500DailyRetS, main="", col="cornflowerblue")
par(mfrow=c(1,1))
```

## Histograms of daily and monthly returns using the same bins

```{r}
msftHist = hist(msftMonthlyRetS, plot=FALSE, breaks=15)
class(msftHist)
names(msftHist)
par(mfrow=c(2,2))
  hist(msftMonthlyRetS, main="", col="cornflowerblue")
  hist(msftDailyRetS, main="", col="cornflowerblue",
       breaks=msftHist$breaks)
  hist(sp500MonthlyRetS, main="", col="cornflowerblue", 
       breaks=msftHist$breaks)
  hist(sp500DailyRetS, main="", col="cornflowerblue",
       breaks=msftHist$breaks)
par(mfrow=c(1,1))

```

## Are Msft returns normally distributed? A first look

```{r}
# create simulated iid Gaussian data with same mean and SD
# as MSFT for both monthly and daily simple returns
set.seed(123)
gwnDaily = rnorm(length(msftDailyRetS), mean=mean(msftDailyRetS),
                 sd=sd(msftDailyRetS))
gwnDaily = xts(gwnDaily, index(msftDailyRetS))
colnames(gwnDaily) = "GWN"
gwnMonthly = rnorm(length(msftMonthlyRetS), mean=mean(msftMonthlyRetS),
                   sd=sd(msftMonthlyRetS))
gwnMonthly = xts(gwnMonthly, index(msftMonthlyRetS))
colnames(gwnMonthly) = "GWN"
```

## Msft monthly returns vs. normal distribution

```{r, echo=FALSE}
# shows time plots and histograms together for Microsoft monthly returns and GWN
par(mfrow=c(2,2))
  plot.zoo(msftMonthlyRetS, main="Monthly Returns on MSFT", 
       lwd=2, col="blue", ylim=c(-0.4, 0.4))
  abline(h=0)
  plot.zoo(gwnMonthly, main="Simulated Normal Returns",  
       lwd=2, col="blue", ylim=c(-0.4, 0.4))
  abline(h=0)
  hist(msftMonthlyRetS, main="", col="cornflowerblue", 
       xlab="returns")
  hist(gwnMonthly, main="", col="cornflowerblue", 
       xlab="returns", breaks=msftHist$breaks)
par(mfrow=c(1,1))
```

## Msft daily returns vs. normal distribution

```{r, echo=FALSE}
msftDailyHist = hist(msftDailyRetS, plot=FALSE, breaks=15)
par(mfrow=c(2,2))
  plot.zoo(msftDailyRetS, main="Monthly Returns on MSFT", 
       lwd=2, col="blue", ylim=c(-0.15, 0.15))
  abline(h=0)
  plot.zoo(gwnDaily, main="Simulated Normal Returns",  
       lwd=2, col="blue", ylim=c(-0.15, 0.15))
  abline(h=0)
  hist(msftDailyRetS, main="", col="cornflowerblue", 
       xlab="returns")
  hist(gwnDaily, main="", col="cornflowerblue", 
       xlab="returns", breaks=msftDailyHist$breaks)
par(mfrow=c(1,1))
```


## Smoothed histogram for Msft monthly returns

```{r}
MSFT.density = density(msftMonthlyRetS)
hist(msftMonthlyRetS, main="", xlab="Microsoft Monthly Returns", 
     col="cornflowerblue", probability=TRUE, ylim=c(0,5.5))
points(MSFT.density,type="l", col="orange", lwd=2)
```

## Empirical CDF of Msft and Sp500 monthly returns

Compute empirical CDF of MSFt returns at a single point

```{r}
# empirical CDF of MSFT returns at R = 0
Fhat.0 = sum(msftMonthlyRetS <= 0)/nrow(msftMonthlyRetS)
Fhat.0
```


```{r}
# compute and plot empirical CDF for all data points
plot(ecdf(coredata(msftMonthlyRetS)), main="", col="blue")

```


## Empirical quantiles of Msft and Sp500 monthly returns

```{r}
# empirical quantiles for MSFT and Sp500
quantile(msftMonthlyRetS)
quantile(sp500MonthlyRetS)
# 1% and 5% quantiles
quantile(msftMonthlyRetS,probs=c(0.01,0.05))
quantile(sp500MonthlyRetS,probs=c(0.01,0.05))
# median and IQR
apply(msftSp500MonthlyRetS, 2, median)
apply(msftSp500MonthlyRetS, 2, IQR)
```

## Using empirical quantiles to compute historical VaR

```{r}
W = 100000
msftQuantiles = quantile(msftMonthlyRetS, probs=c(0.01, 0.05))
sp500Quantiles = quantile(sp500MonthlyRetS, probs=c(0.01, 0.05))
msftVaR = W*msftQuantiles
sp500VaR = W*sp500Quantiles
msftVaR
sp500VaR
```

## Normal QQ-plots for GWN, Microsoft and S&P 500 returns

```{r}
par(mfrow=c(2,3))  # 4 panel layout: 2 rows and 2 columns
	qqnorm(gwnMonthly, main="GWN Monthly", col="cornflowerblue")
	qqline(gwnMonthly)
	qqnorm(msftMonthlyRetS, main="MSFT Monthly Returns", col="cornflowerblue")
	qqline(msftMonthlyRetS)
	qqnorm(sp500MonthlyRetS, main="SP500 Monthly Returns", col="cornflowerblue")
	qqline(sp500MonthlyRetS)
  qqnorm(gwnDaily, main="GWN Daily", col="cornflowerblue")
  qqline(gwnDaily)
  qqnorm(msftDailyRetS, main="MSFT Daily Returns", col="cornflowerblue")
  qqline(msftDailyRetS)
  qqnorm(sp500DailyRetS, main="SP500 Daily Returns", col="cornflowerblue")
  qqline(sp500DailyRetS)
par(mfrow=c(1,1))
```

## Student's t QQ-plot for Msft returns

```{r}
library(car)
qqPlot(coredata(msftMonthlyRetS), distribution="t", df=5, 
       ylab="MSFT quantiles", envelope=FALSE)
```

## Sample shape statistics for the monthly returns on Msft and Sp500

```{r}
# Monthly return stats
statsMonthly = rbind(apply(msftSp500MonthlyRetS, 2, mean),
                     apply(msftSp500MonthlyRetS, 2, var),
                     apply(msftSp500MonthlyRetS, 2, sd),
                     apply(msftSp500MonthlyRetS, 2, skewness),
                     apply(msftSp500MonthlyRetS, 2, kurtosis))
rownames(statsMonthly) = c("Mean", "Variance", "Std Dev",
                           "Skewness",  "Excess Kurtosis")
round(statsMonthly, digits=4)
```

The annualized mean and volatility are

```{r}
12*statsMonthly["Mean", ]
sqrt(12)*statsMonthly["Std Dev", ]
```


## Sample shape statistics for the monthly returns on Msft and Sp500

```{r}
# Daily return stats
statsDaily = rbind(apply(msftSp500DailyRetS, 2, mean),
                   apply(msftSp500DailyRetS, 2, var),
                   apply(msftSp500DailyRetS, 2, sd),
                   apply(msftSp500DailyRetS, 2, skewness),
                   apply(msftSp500DailyRetS, 2, kurtosis))
rownames(statsDaily) = rownames(statsMonthly)
round(statsDaily, digits=4)
```

Check the square root of time rule (20 trading days in a month)

```{r}
20*statsDaily["Mean", ]
sqrt(20)*statsDaily["Std Dev", ]
```

## Outliers

```{r, echo=FALSE}
# create nice systematic relationship with an outlier
n.obs = 20
set.seed(123)
x = seq(0, 10, length.out=n.obs)
e = rnorm(n.obs, 0, 0.1)
y = 1 + 0.5*x + e
y[5] = 5
plot(x, y, type="p", lwd=2, col="blue", pch=16, cex=1.5)
points(x[5], y[5], pch=16, cex=1.5)
```


## Effect of outliers on sample statistics

```{r}
# create GWN return data polluted by outlier
gwnMonthlyOutlier = gwnMonthly
gwnMonthlyOutlier[20] = -sd(gwnMonthly)*6

par(mfrow=c(2,1))
  plot.zoo(gwnMonthlyOutlier,main="", lwd=2, col="blue", ylab="")
	abline(h=0)
	hist(gwnMonthlyOutlier, main="", col="cornflowerblue", 
       xlab="GWN with outlier")
par(mfrow=c(1,1))

gwnMonthlyBoth = cbind(gwnMonthly, gwnMonthlyOutlier)
statsMat = rbind(apply(gwnMonthlyBoth, 2, mean),
                 apply(gwnMonthlyBoth, 2, var),
                 apply(gwnMonthlyBoth, 2, sd),
                 apply(gwnMonthlyBoth, 2, skewness),
                 apply(gwnMonthlyBoth, 2, kurtosis),
                 apply(gwnMonthlyBoth, 2, median),
                 apply(gwnMonthlyBoth, 2, IQR))
rownames(statsMat) = c("Mean", "Var", "SD", "skewness", "kurtosis",
                       "median", "IQR")
statsMat["kurtosis", ] = statsMat["kurtosis", ] + 3
pctchange = (statsMat[, 2] - statsMat[, 1])/statsMat[, 1]
statsMat = cbind(statsMat, pctchange)
round(statsMat, digits=4)
```

## Boxplots

A Boxplot gives a graphical summary of a data distribution using outlier robust statistics

```{r, echo=FALSE}
q.25 = quantile(msftMonthlyRetS, 0.25)
q.75 = quantile(msftMonthlyRetS, 0.75)
q.5 = quantile(msftMonthlyRetS, 0.5)
iqr = IQR(msftMonthlyRetS)
boxplot(coredata(msftMonthlyRetS), col="cornflowerblue")
text(x=1.25, y = q.25, labels="q.25")
text(x=1.25, y = q.5, labels="q.50")
text(x=1.25, y = q.75, labels="q.75")
text(x=1.25, y = q.25-iqr, labels="q.25 - 1.5*IQR")
text(x=1.25, y = q.75+iqr, labels="q.75 + 1.5*IQR")
text(x=1.25, y = 0.3, labels="extreme values")
text(x=1.25, y = -0.3, labels="extreme values")
```


## Boxplots of return distributions

Boxplots are very good for comparing the distributions of multiple asset returns:

```{r}
dataToPlot = merge(msftMonthlyRetS, msftMonthlyRetC, sp500MonthlyRetS, sp500MonthlyRetC)
colnames(dataToPlot) = c("msftRetS", "msftRetC", "sp500RetS", "sp500RetC")
boxplot(coredata(dataToPlot), main="Simple vs. CC Monthly Returns",
        col="cornflowerblue")
```        

## Time series descriptive statistics

Sample autocorrelations

```{r}
acf(coredata(msftMonthlyRetS), lag.max=5, plot=FALSE)
acf(coredata(msftDailyRetS), lag.max=5, plot=FALSE)

```

## SACF

```{r}
par(mfrow=c(2,2))
  acf(coredata(msftMonthlyRetS), main="msftMonthlyRetS", lwd=2)
  acf(coredata(sp500MonthlyRetS), main="sp500MonthlyRetS", lwd=2)
  acf(coredata(msftDailyRetS), main="msftDailyRetS", lwd=2)
  acf(coredata(sp500DailyRetS), main="sp500DailyRetS", lwd=2)
par(mfrow=c(1,1))
```

## Time dependence in volatility

```{r}
n.obs = 500
set.seed(123)
gwn = rnorm(n.obs)
```

```{r, echo=FALSE}
par(mfrow=c(3, 2))
  ts.plot(gwn)
  abline(h=0)
  acf(gwn, lwd=2)
  ts.plot(gwn^2)
  acf(gwn^2, lwd=2)
  ts.plot(abs(gwn))
  acf(abs(gwn), lwd=2)
par(mfrow=c(1,1))
```

## Time dependence in S&P 500 daily volatility

```{r}
dataToPlot = merge(sp500DailyRetS, abs(sp500DailyRetS), sp500DailyRetS^2)
colnames(dataToPlot) = c("Returns", "abs(Returns)", "Returns^2")
plot.zoo(dataToPlot, main="Daily Returns", col="blue")
```

```{r}
par(mfrow=c(2,2))
acf(abs(coredata(sp500DailyRetS)), main="Daily Absolute Returns", lwd=2)
acf(coredata(sp500DailyRetS)^2, main="Daily Squared Returns", lwd=2)
acf(abs(coredata(sp500MonthlyRetS)), main="Monthly Absolute Returns", lwd=2)
acf(coredata(sp500MonthlyRetS)^2, main="Monthly Squared Returns", lwd=2)
par(mfrow=c(1,1))
```

## Bivariate descriptive statistics

```{r}
par(mfrow=c(1,2))
plot(coredata(sp500MonthlyRetS),coredata(msftMonthlyRetS), 
         main="Monthly returns", 
         xlab="S&P500", ylab="MSFT",
         lwd=2, pch=16, cex=1.25, col="blue")
abline(v=mean(sp500MonthlyRetS))  
abline(h=mean(msftMonthlyRetS))
plot(coredata(sp500DailyRetS),coredata(msftDailyRetS), 
         main="Daily returns", 
         xlab="S&P500", ylab="MSFT",
         lwd=2, pch=16, cex=1.25, col="blue")
abline(v=mean(sp500DailyRetS))  
abline(h=mean(msftDailyRetS))
```

## Pairs plot

```{r}

pairs(coredata(merge(gwnMonthly,msftMonthlyRetS,sp500MonthlyRetS)), col="blue", 
      pch=16, cex=1.25, cex.axis=1.25)
```

## Sample covariance and correlations

For Monthly data, the sample covariance and correlation between Microsoft and the S&P 500 simple returns are:

```{r}
cov(sp500MonthlyRetS, msftMonthlyRetS)
cor(sp500MonthlyRetS, msftMonthlyRetS)
```

The sample covariance and correlation matrices are

```{r}
cov(msftSp500MonthlyRetS)
cor(msftSp500MonthlyRetS)
```


For daily data, the sample covariance and correlation between Microsoft and the S&P 500 simple returns are:

```{r}
cov(sp500DailyRetS, msftDailyRetS)
cor(sp500DailyRetS, msftDailyRetS)
```

The sample covariance and correlation matrices are

```{r}
cov(msftSp500DailyRetS)
cor(msftSp500DailyRetS)
```

Check aggregation (square-root-of-time) rule used for variances

```{r}
20*cov(sp500DailyRetS, msftDailyRetS)
```

## Sample cross-lag covariances and correlations

```{r}
Ghat = acf(coredata(msftSp500MonthlyRetS), type="covariance", lag.max=5, plot=FALSE)
Chat = acf(coredata(msftSp500MonthlyRetS), type="correlation", lag.max=5, plot=FALSE)
class(Ghat)
names(Ghat)
```

Extract the cross-lag correlation matrices C0 and C1:

```{r}
# Chat0
Chat$acf[1,,]
# Chat1
Chat$acf[2,,]
```

Show sample autocorrelations and all pairwise cross-lag correlations:

```{r}
Chat
```

Plot these values

```{r}
plot(Chat, lwd=2)
```

