print(subtime[ind])
print(x.fit[[1]])
print(x.fit[[2]])
time<-rbind(time,as.character(subtime[ind]))
fitin<-rbind(fitin,c(as.matrix(x.fit[[1]])[2]+as.matrix(x.fit[[1]])[3],as.matrix(x.fit[[1]])[2],as.matrix(x.fit[[1]])[3]))
fitout<-rbind(fitout,c(as.matrix(x.fit[[2]])[2]+as.matrix(x.fit[[2]])[3],as.matrix(x.fit[[2]])[2],as.matrix(x.fit[[2]])[3]))
indexnew = ((1+(i-1)*1000)+ ind)
#print(timestamp)
indexold = c(indexold, indexnew)
}
}
}
###Results in Table 2 will be printed afterwards
###Store results in latex
result_cz<-data.frame(time,fitin,fitout,sta)
result_cz<-result_cz[2:6,]
rownames(result_cz)<-c(1:nrow(result_cz))
colnames(result_cz)<-c("time1","time2","in1","in2","in3","out1","out2","out3","sta")
result_cz
print.xtable(xtable(head(result_cz)), file = "./Table/table2.tex")
indices = as.vector(as.numeric(indexold))
indexvix = indices[-1]
# return the indices of change point. this is saved for the use of movingwindowfit_vix_2
save(indexvix, file = "indexvix.RData")
################################################################
#                                                              #
#    This Program Plots Figure 9 in the Paper                  #
#   You will need to run "garch_vix_new.r" first               #
#                                                              #
################################################################
rm(list=ls())
###Setting working directory to where you store the replication package folder
#e.g. setwd("C:/Users/Desktop/submission_final/3Replication_Package")
setwd(changing to the path where you store the replication package folder)
library("aTSA")
library("fBasics")
library("truncnorm")
library("rugarch")
library("zoo")
library("lmtest")
library("forecast")
library("xts")
library("tseries")
library("fGarch")
library("TTR")
source("helpfun.r")
library("Rcpp")
Rcpp::sourceCpp('garchlikelihood.cpp')
Rcpp::sourceCpp('garchlikelihoodgradient.cpp')
###Read the data "indexvix.RData", that contains the estimation results by running "garch_vix_new".
load("indexvix.RData")
### Also the TED data, named "Vix.txt" in the replication package folder.
Vix = read.table("Data/Vix.txt", header = FALSE, sep = "", dec = ".")
Vix = Vix[-1,]
time.stamp = Vix[,1]
nn          = nrow(Vix)
Vixm       = matrix(as.numeric(unlist(Vix[,2:5])),nn,4)
Vixm       = Vixm[,4]
dvix        = diff(log(Vixm))
adf.test(dvix,output = TRUE)
library("aTSA")
library("fBasics")
library("truncnorm")
library("rugarch")
library("zoo")
library("lmtest")
library("forecast")
library("xts")
library("tseries")
library("fGarch")
library("TTR")
source("helpfun.r")
library("Rcpp")
Rcpp::sourceCpp('garchlikelihood.cpp')
Rcpp::sourceCpp('garchlikelihoodgradient.cpp')
###Read the data "indexvix.RData", that contains the estimation results by running "garch_vix_new".
load("indexvix.RData")
### Also the TED data, named "Vix.txt" in the replication package folder.
Vix = read.table("Data/Vix.txt", header = FALSE, sep = "", dec = ".")
Vix = Vix[-1,]
time.stamp = Vix[,1]
nn          = nrow(Vix)
Vixm       = matrix(as.numeric(unlist(Vix[,2:5])),nn,4)
Vixm       = Vixm[,4]
dvix        = diff(log(Vixm))
###adf.test(dvix,output = TRUE)
###pacf(dvix)
###acf(dvix)
fit = arima(dvix, order = c(9,0,2))
###pacf(fit$residuals)
xx = dvix
plot(xx^2,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
timestamp = as.Date(strptime(time.stamp, "%m/%d/%Y"))
x = xx
index = indexvix
###index = c(562,802,2,282,1702,1822,1121,1602,1862,1982,2002,2282)
ll    = length(index)/2
indexsort = sort(index)
indexeven  = seq(1,ll)*2
indexodd   = seq(1,ll)*2-1
basicStats(x)
Sys.setlocale("LC_TIME", "C")
pdf("vix.pdf", width=12, height=6)
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
plot(timestamp,Vixm,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"), format="%B %Y") )
dev.off()
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
plot(timestamp,Vixm,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"), format="%B %Y") )
library("aTSA")
library("fBasics")
library("truncnorm")
library("rugarch")
library("zoo")
library("lmtest")
library("forecast")
library("xts")
library("tseries")
library("fGarch")
library("TTR")
source("helpfun.r")
library("Rcpp")
Rcpp::sourceCpp('garchlikelihood.cpp')
Rcpp::sourceCpp('garchlikelihoodgradient.cpp')
###Read the data "index.RData", that contains the estimation results by running "garch_bitcoin_app".
load("index.RData")
#read the data
bit  = read.table("Data/bitcoin.txt", header = FALSE, sep = "", dec = ".")
bit = bit[-1,]
time.stamp = bit[,1]
timestamp = as.Date(strptime(time.stamp, "%d/%m/%Y"))
nn          = nrow(bit)
Vixm       = matrix(as.numeric(unlist(bit[,3])),nn,1)
dvix        = diff(log(Vixm))
pacf(dvix)
acf(dvix)
adf.test(dvix)
fit = arima(dvix, order = c(6,0,0))
pacf(fit$residuals)
x = dvix
pdf("bit.pdf", width=12, height=6)
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
timestamp = as.Date(strptime(time.stamp, "%d/%m/%Y"))
plot(timestamp,Vixm,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"), format="%B %Y") )
dev.off()
#plot the data with the change point
ll    = length(index)/2
indexsort = sort(index)
indexeven  = seq(1,ll)*2
indexodd   = seq(1,ll)*2-1
basicStats(x)
Sys.setlocale("LC_TIME", "C")
#moving window fit of the parameters
x = fit$residuals*20
kernliknuold2 <- function(a,x,r1,r2) {
return(garchlikelihood(a,x,r1,r2));
}
kernliknuold2grad <- function(a,x,r1,r2) {
return(garchlikelihoodgradient(a,x,r1,r2));
}
kernliknuold2aussen <- function(a,x,r1,r2) {
return(garchlikelihood(a,x,1,r1) + garchlikelihood(a,x,r2,n));
}
kernliknuold2gradaussen <- function(a,x,r1,r2) {
return(garchlikelihoodgradient(a,x,1,r1) + garchlikelihoodgradient(a,x,r2,n));
}
#x = x[1:1000]
n = length(x)
windowsize = 200
endwin     = length(x) - windowsize
paramatrix = matrix(0,(endwin),3)
for (i in 1:endwin)
{
x.fit  = estimate(c(i,i+windowsize),x)
paramatrix[i,] = x.fit[[1]]
print(x.fit[[1]])
}
abline(v=15, col="blue")
###Figure 1 in the paper
pdf("Figures/paravarying_bit.pdf", width=12, height=6)
thetan     =  paramatrix
xx = seq(1, length(thetan[,2]))
xx = seq(1:i)
#xx  =  seq((windowsize+1):length(x))
xxt = timestamp[xx]
plot(xxt,thetan[,2]+ thetan[,3],type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "", ylim = c(0,1.5),lty = 6)
lines(xxt,thetan[,2],type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
lines(xxt,thetan[,3],type= "l",col = "black",xaxt="n", lwd = 2, xlab = "",ylab=  "",lty = 3)
abline(h= 1, col="grey")
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"),format="%B %Y"),cex.axis=1 )
dev.off()
library("aTSA")
library("fBasics")
library("truncnorm")
library("rugarch")
library("zoo")
library("lmtest")
library("forecast")
library("xts")
library("tseries")
library("fGarch")
library("TTR")
source("helpfun.r")
library("Rcpp")
library("xtable")
Rcpp::sourceCpp('garchlikelihood.cpp')
Rcpp::sourceCpp('garchlikelihoodgradient.cpp')
### Reading the VIX data, named "Vix.txt" in the replication package folder.
indexold = c(0)
Vix = read.table("Data/Vix.txt", header = FALSE, sep = "", dec = ".")
Vix = Vix[-1,]
time.stamp = Vix[,1]
nn          = nrow(Vix)
Vixm       = matrix(as.numeric(unlist(Vix[,2:5])),nn,4)
Vixm[,4]
timestamp = as.Date(strptime(time.stamp, "%m/%d/%Y"))
plot(timestamp,Vixm[,4],type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
lines(timestamp,Vixm[,4],lwd = 2,col = "red")
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"),"%y"),las=2 )
Vixm       = Vixm[,4]
dvix        = diff(log(Vixm))
adf.test(dvix)
pacf(dvix)
acf(dvix)
fit = arima(dvix, order = c(9,0,2))
pacf(fit$residuals)
xx = dvix
plot(xx^2,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
timestamp = as.Date(strptime(time.stamp, "%m/%d/%Y"))
x = xx
basicStats(x)
### Plot of Figure 7 in the paper
pdf("Figures/vix.pdf", width=12, height=6)
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
plot(timestamp,Vixm,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"),"%y"),las=2,cex.axis=1 )
dev.off()
####Plot of Figure 8 in the paper
pdf("Figures/qqvix.pdf", width=12, height=6)
layout(matrix(c(1,2),1,2, byrow = TRUE))
par(mar = c(2,2,1,1) )
hist(x, xlab="Daily return of VIX index", prob=TRUE, main="Histogram for daily return of the VIX index", ylim = c(0,6),col = NULL)
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
lines(xfit, yfit, col="blue", lwd=2)
qqnorm(x)
qqline(x, col = 2)
dev.off()
#yy = fit$residuals
yy = fit$residuals
yy = x*20
###To store the estimation results
sta<-0
time<-0
fitin<-0
fitout<-0
kernliknuold2 <- function(a,x,r1,r2) {
return(garchlikelihood(a,x,r1,r2));
}
kernliknuold2grad <- function(a,x,r1,r2) {
return(garchlikelihoodgradient(a,x,r1,r2));
}
kernliknuold2aussen <- function(a,x,r1,r2) {
return(garchlikelihood(a,x,1,r1) + garchlikelihood(a,x,r2,n));
}
kernliknuold2gradaussen <- function(a,x,r1,r2) {
return(garchlikelihoodgradient(a,x,1,r1) + garchlikelihoodgradient(a,x,r2,n));
}
win = 1000
result= list()
#yy = x
for (i in 1: floor(length(yy)/win))
{
# i = 2
if (i ==3)
{x = yy[(1+(i-1)*1000):length(yy)]
subtime = timestamp[(1+(i-1)*1000):length(timestamp)]
}
if (i!=3)
{
x = yy[(1+(i-1)*1000):(i*1000)]
subtime = timestamp[(1+(i-1)*1000):(i*1000)]
}
n = length(x)
x = scalex(x,c(180,90))
###
grid = 20
kappa=0.1;
kappastrich=0.1;
chi = 0.5;
H    = c(0,1,1)
#Running index tau1
indizes_tau1 = seq(1,n,by=grid);
indizes_tau2 = seq(1,n,by=grid);
be = matrix(0,nrow=n,ncol=n);
x.fit2 = garchFit(~ garch(1,1), data = x, trace = FALSE, cond = "QMLE",include.mean= FALSE)
rang = 800
for (tau1n in indizes_tau1 )
{
indizes_tau2 = indizes_tau1[indizes_tau1 - tau1n > floor(kappa*n) & floor((1-kappastrich)*n) > indizes_tau1 - tau1n & indizes_tau1 - tau1n < rang];
if (length(indizes_tau2) > 0){
for (tau2n in indizes_tau2) {
##  x.fit  = garchFit(~ garch(1,1), data = x[tau1n:tau2n], trace = FALSE, cond = "QMLE",include.mean= FALSE)
x.fit  = estimate(c(tau1n,tau2n),x)
##xnew   = c(x[1:tau1n],x[tau2n:n])
###  tt     = c(1,(length(xnew)-1))
###x.fit2 = estimate(tt,xnew)
be[tau1n,tau2n] = ((tau2n-tau1n)/(n))^{chi}*sqrt((solve(t(H)%*%x.fit[[3]]%*%H)))* sqrt(n)*(t(H)%*%x.fit[[1]]-0.95)
# print(t(H)%*%x.fit[[1]])
# if( be[tau1n,tau2n]>3.29)
#{print(t(H)%*%x.fit[[1]])}
}
}
}
tsmax = max((be))
amax = which(be==max(be), arr.ind=T)
CV = 3.2
if(tsmax  > 3.2)
{
stampall = which(be>= CV, arr.ind=T)
stampmax = which(be==max(be), arr.ind=T)
}
benew = be
beneww = be
stampmaxall = c(0,0)
while(length(which(beneww> CV, arr.ind=T))!=0)
{
stampallnew  = which(beneww> CV, arr.ind=T)
stampmaxnew  = which(beneww==max(beneww), arr.ind=T)
print(max(beneww))
sta<-rbind(sta,max(beneww))
#print(stampmaxall)
colindex           = which((stampallnew[,1]>stampmaxnew[2]))
rowindex           = which(stampallnew[,2]< stampmaxnew[1])
indices            = union(colindex,rowindex)
stampallnew        = stampallnew[indices,]
beneww              = matrix(0,n,n)
beneww[stampallnew] = benew[stampallnew]
#print(subtime[stampallnew])
stampmaxall = rbind( stampmaxall, stampmaxnew)
}
print(i)
stampmaxall = stampmaxall[-1,]
print(stampmaxall)
nstamp = nrow((stampmaxall))
if(is.null(ncol(stampmaxall)))
{
nstamp = 1
h0=0.95
ind  = stampmaxall
x.fit  = estimate_new(ind,x,h0)
print(subtime[ind])
print(x.fit[[1]])
print(x.fit[[2]])
time<-rbind(time,as.character(subtime[ind]))
fitin<-rbind(fitin,c(as.matrix(x.fit[[1]])[2]+as.matrix(x.fit[[1]])[3],as.matrix(x.fit[[1]])[2],as.matrix(x.fit[[1]])[3]))
fitout<-rbind(fitout,c(as.matrix(x.fit[[2]])[2]+as.matrix(x.fit[[2]])[3],as.matrix(x.fit[[2]])[2],as.matrix(x.fit[[2]])[3]))
indexnew = ((1+(i-1)*1000)+ ind)
indexold = c(indexold, indexnew)
} else  {
for (ss in 1:nstamp)
{
h0=0.95
ind  = stampmaxall[ss,]
x.fit  = estimate_new(ind,x,h0)
print(subtime[ind])
print(x.fit[[1]])
print(x.fit[[2]])
time<-rbind(time,as.character(subtime[ind]))
fitin<-rbind(fitin,c(as.matrix(x.fit[[1]])[2]+as.matrix(x.fit[[1]])[3],as.matrix(x.fit[[1]])[2],as.matrix(x.fit[[1]])[3]))
fitout<-rbind(fitout,c(as.matrix(x.fit[[2]])[2]+as.matrix(x.fit[[2]])[3],as.matrix(x.fit[[2]])[2],as.matrix(x.fit[[2]])[3]))
indexnew = ((1+(i-1)*1000)+ ind)
#print(timestamp)
indexold = c(indexold, indexnew)
}
}
}
###Results in Table 2 will be printed afterwards
###Store results in latex
result_cz<-data.frame(time,fitin,fitout,sta)
result_cz<-result_cz[2:6,]
rownames(result_cz)<-c(1:nrow(result_cz))
colnames(result_cz)<-c("time1","time2","in1","in2","in3","out1","out2","out3","sta")
result_cz
print.xtable(xtable(head(result_cz)), file = "./Table/table2.tex")
indices = as.vector(as.numeric(indexold))
indexvix = indices[-1]
# return the indices of change point. this is saved for the use of movingwindowfit_vix_2
save(indexvix, file = "indexvix.RData")
###To run the code successfully, you will also first need to install the following packages
###These packages only need to be installed once. and you may ignore this if they are already installed in your R/Rstudio.
#install.packages("aTSA","fBasics","truncnorm","rugarch","zoo","lmtest","forecast","xts","tseries","fGarch","TTR","Rcpp")
#######################################################################################
#                                                                                     #
#    Changing path and installing packages are the two only things you need to do     #
#    to run this programme after downloading the whole replication package            #
#                                                                                     #
#######################################################################################
###Load all the following packages as follows and also source some help functions.
###The files "helpfun.r", "garchlikelihood.cpp","garchlikelihoodgradient.cpp" contain functions used in this code.
#These files are provided in the replication package folder.
library("aTSA")
library("fBasics")
library("truncnorm")
library("rugarch")
library("zoo")
library("lmtest")
library("forecast")
library("xts")
library("tseries")
library("fGarch")
library("TTR")
source("helpfun.r")
library("Rcpp")
Rcpp::sourceCpp('garchlikelihood.cpp')
Rcpp::sourceCpp('garchlikelihoodgradient.cpp')
###Read the data "indexted.RData", that contains the estimation results by running "garch_ted_spread_2021_02_17".
load("indexted.RData")
index = indexred
ll    = length(index)/2
indexsort = sort(index)
indexeven  = seq(1,ll)*2
indexodd   = seq(1,ll)*2-1
Sys.setlocale("LC_TIME", "C")
### Also the TED data, named "ted_spread.txt" in the replication package folder.
Vix        = read.table("Data/ted_spread.txt", header = FALSE, sep = "", dec = ".")
time.stamp = Vix[,1]
nn         = nrow(Vix)
Vixm       = matrix(as.numeric(unlist(Vix[,2])),nn,1)
Vixm
time.stamp.new =  time.stamp[Vixm != 0]
length(time.stamp)
Vixm_new       = Vixm[Vixm != 0]
timestamp = as.Date(strptime(time.stamp.new, "%Y-%m-%d"))
plot(timestamp,Vixm_new ,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
lines(timestamp,Vixm_new ,lwd = 2,col = "red")
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"),"%y"),las=2 )
dvix        = diff(log(Vixm_new))
x = dvix
###
pdf("ted.pdf", width=12, height=6)
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
plot(timestamp,Vixm_new,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"), format="%B %Y") )
dev.off()
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
plot(timestamp,Vixm_new,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"), format="%B %Y") )
pdf("ted.pdf", width=12, height=6)
layout(matrix(c(1,2,3),3,1, byrow = TRUE))
par(mar = c(2,2,1,1) )
plot(timestamp,Vixm_new,type= "l",col = "blue",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],x,type= "l",col = "red",xaxt="n", lwd = 2, xlab = "",ylab=  "")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
plot(timestamp[-1],abs(x),type='l', lwd = 2, xaxt="n",xlab = "",ylab=  "", main="")
abline(v = timestamp[index[indexeven]], col = "lightgray", lty = 6,lwd = 4)
abline(v = timestamp[index[indexodd]], col = "darkgray", lty = 6,lwd = 4)
axis.Date(1,at=seq(min(timestamp), max(timestamp), by="year"),labels=format(seq(min(timestamp), max(timestamp), by="year"), format="%B %Y") )
dev.off()
