# R 3.4.1
# daily, monthly returns: plot, distribution Ex3_month_daily_VRtest_R.txt
# VR test for white noise
# resampling distribution for VR test

library("tseries")
library("zoo")
library("moments")
 
####################################
###  script for basic statistics ###
####################################
basic_stats <- function(ts0) {
     ts <- as.numeric(ts0)
     hh <- c("Basic statistics:    "," Nobs =", formatC(length(ts), width=9) )
     print.noquote(hh)
     hh <- c("Mean=    ", formatC( mean(ts), digits = 5, width= 8, format = "fg") , 
             "StDev=   ", formatC(  sd(ts) , digits = 5, width= 8, format = "fg")  )
     print.noquote(hh )                  
     hh <- c("Skewness=",formatC( skewness(ts), digits = 5, width= 8, format = "fg") ,
             "Kurtosis=",formatC( kurtosis(ts), digits = 5, width= 8, format = "fg")  )
     print.noquote(hh) 
     jarque.bera.test(ts)
                           }  
#####################################

start0 <- "1998-10-01"
end0   <- "2004-10-01"  

#######################################
### compare two indices, daily data ###
#######################################
# Read data
require(quantmod) # Finanzdaten download

spc <- new.env()
dax <- new.env()

setDefaults(getSymbols,src="yahoo")
getSymbols("^GSPC",  env = spc, from = start0, to = end0 )
getSymbols("^GDAXI", env = dax, from = start0, to = end0 )
# read in series are of xts type,  plot is different for xts and ts

x.d <-as.zoo( spc$GSPC[,4]);       y.d <- as.zoo(dax$GDAXI[,4]); 

## Plot series, returns 
plot(merge(x.d,y.d), main = "Index: ^GSPC, ^GDAXI, daily")
#plot(x.d, y.d) 
r.xd  <- diff(log(x.d));         # return US stock index
r.yd  <- diff(log(y.d));         # return other index
plot( merge(r.xd, r.yd))         # when US is coughing, EU has fever

r.xd <- r.xd[is.na(r.xd)==FALSE]
r.yd <- r.yd[is.na(r.yd)==FALSE]

## Distribution of returns, 
plot(density(r.xd), main="r-^GSPC ... black, r-^GDAXI ... blue")
lines(density(r.yd), col="blue")  # 2nd plot printed in 1st 

# basic statistics
basic_stats(r.xd)
basic_stats(r.yd)

#########################################
### compare two indices, monthly data ###
#########################################
### monthly data are the end-of-month value:  index is in days!
x.m <- x.d[!duplicated(as.yearmon(time(x.d)), fromLast = TRUE)]
y.m <- y.d[!duplicated(as.yearmon(time(y.d)), fromLast = TRUE)]
                                                 
#plot(x.m);  plot(y.m);                  
                            # merge works on daily bases 
                            # end-of-month-day different in both indices
xm <- ts(as.numeric(x.m), start = c(1998, 10), frequency = 12 )
ym <- ts(as.numeric(y.m), start = c(1998, 10), frequency = 12 )
                            # merge does not work for ts
x.mz <- as.zoo(xm)
y.mz <- as.zoo(ym)
plot(merge(x.mz,y.mz), main = "Index: ^GSPC, ^GDAXI, monthly")
#plot(x.mz, y.mz) 

r.xmz  <- diff(log(x.mz));         # return US stock index
r.ymz  <- diff(log(y.mz));         # return other index

#plot(r.xmz); plot(r.ymz);
plot( merge(r.xmz, r.ymz), main = "Return: ^GSPC, ^GDAXI, monthly")      

r.xmz <- r.xd[is.na(r.xmz)==FALSE]
r.ymz <- r.yd[is.na(r.ymz)==FALSE]

## Distribution of returns, 
plot(density(r.xmz), main="r-^GSPC ... black, r-^GDAXI ... blue")
lines(density(r.ymz), col="blue")        # 2nd plot printed in 1st 

# basic statistics
basic_stats(r.xmz)
basic_stats(r.ymz)

###########################
### Variance Ratio test ###    # based only on daily data 
###########################    # a month is assumed to have 22 trading days
p <- as.vector(log(x.d));       

m <- 22;
M <- as.integer( (length(p)-1) / m ) 
TT <- m*M
p <- p[1:(TT+1)];        # truncate series to full months

## calculation of the sample statistic
r <- diff(p);
#r <- rnorm(TT)          # for checking whether resampling yields the asy distr
mu <- mean(r)

rm <- c(1:M)*0
for (i in 1:M) { j2 <- m*i; j1 <- m*(i-1)+1;  sumr <- 0;
                 for (ii in j1:j2) {  sumr <- sumr + r[ii]; }
                 rm[i] <- sumr }      # monthly returns rm

s2a <- var(r)*(TT-1)/TT               # var() ... variance corrected for df's 
mcrm <- rm - m*mu                     # mean corrected monthly returns mcrm
s2b <- ( var(mcrm)*(M-1)/M )/m
VRm <- s2b/s2a -1           # asy distributed as N(0,2*(m-1)/TT) under H0
VRm                        
1.96*sqrt(2*(m-1)/TT)

## resampling of r  ( mu = mean(r), s2a are invariant )
## small sample distribution of VR is simulated
nsim <- 1000
o <- c(1:TT)
VRms <- c(1:nsim)*0

set.seed(1234)
for (is in 1:nsim) 
 {
   o   <- sample(o) 
   rs  <- r[o]  
   rms <- c(1:M)*0
   for (i in 1:M) { j2 <- m*i; j1 <- m*(i-1)+1; sumr <- 0; 
                    for (ii in j1:j2) {  sumr <- sumr + rs[ii]; }
                    rms[i] <- sumr }        # monthly returns rm

    mcrms <- rms - m*mu                     # mean corrected monthly returns mcrm
    s2bs <- ( var(mcrms)*(M-1)/M )/m
    VRms[is] <- s2bs/s2a -1 
 }

VRm
quantile(VRms, c(0.05, 0.95) ) 

## comparison, asy and resampling distribution
plot(function(x) dnorm( x,0,(2*(m-1)/TT)^(0.5) ), -1, 1, main="empirical density ... blue ") 
lines(density(VRms), col="blue") 
     
## EXERCISE:
#  (1) Read in daily and monthly returns and compare their distrubutions.
#  (2) Test for white noise with the VR test.
#  (3) Calculate the empirical distribution of the VR test via the resmapling method 
#      and use it for testing.