# R 3.0.1
# distribution of intraday/realized volatility       Ex3_RealVola_msft_R.txt
#
# read msft 1-minute tick data 08/05/2008 to 26/11/2008 (dd/mm/yyyy)

#setwd("C:/Users/SmartBoard/Downloads/")  

P <- read.table("Tick_MSFT_V1_1c")      # change working directory
colnames(P) <- c("day","min","price","dm","hhm","d")
head(P)
tail(P)
 # day, min ... date, dm ... minute in a single day (1-391) 
 # hhm ... 1/2-hour minute (0,1-30), d ... day (1-140)
 # trading hours NY 9:30 - 16:00 == 13 (1/2) hours == 391 minutes per day
 # (time stamp is local time)
 # 140 days in total
 # opening price = (dm=1/hhm=0)

# nobs = 54740 = 140 * [ 13      * 30        + 1 ]
#                days  1/2-hours   minutes     Open price
#                                 per 1/2-hour 
# P ... 391 quotes per day, 
# r ... 390 quotes per day + 1 overnight return
Npq <- 391
Nrq <- 390

p = as.ts(log(P[,3]))  # freq =391 does not work; ts plots are nicer

######################
### log price data ###
######################
dip_a <- function(d) { (d-1)*Npq + 1 }          # begin of day / price
dip_e <- function(d) {  dip_a(d) + Npq - 1 }    # end of day  / price
dip_o <- function(d) {  dip_a(d) + Npq }  # overnight price, end of day
                                          #     not available for d=140

pd <- function(d) { p[dip_a(d):dip_e(d)] }      # picks quotes of day d
#pd <- function(d) { p[dip_a(d):dip_o(d)] }      # with overnight price
      
pp <- as.ts( cbind( pd(1), pd(70), pd(140)) ) # for 3 days: 1st, 70th, last
plot(pp)

###################
### log returns ###
###################
r <- diff(p)                                    # 1-minute returns

dir_a <- function(d) { (d-1)*Nrq + (d-1) + 1 }  # begin of day / return
dir_e <- function(d) {  dir_a(d) + Nrq -1  }    # end of day  / return
dir_o <- function(d) {  dir_a(d) + Nrq }  # overnight return, end of day
                                          #     not available for d=140
rd <- function(d) { r[dir_a(d):dir_e(d)] }    # picks returns of day d
#rd <- function(d) { r[dir_a(d):dir_o(d)] }    # with overnight return 
     
rr <- as.ts( cbind( rd(1), rd(70), rd(139)) ) # for 3 days: 1st, 70th, (last-1)
plot(rr)

plot( density(rd(1)), main="within day returns" )
lines( density(rd(70)) , col="blue" )
lines( density(rd(139)) , col="red" )         # compares 3 return distributions

acf(rd(1))                                    # ACF of 1-minute returns of day 1

############################
### rv (within day only) ### 
############################
#      between/overnight only (r^2_o)

mrv <- 5    # choose mrv = 1,2,3,4,5,6,10,12,15,20,30 minutes
Nrvq <- 390/mrv        # no of intervals for rv per (within)day

lrv <- Nrvq*140 + (140-1) # length of rv[.]
rv  <- c(1:lrv)*0

#######
tr <- 0;   trv <- 0;   # time count for r[.] and rv[.] set 0
for (id in 1:140)        
 {                       # for each day
   meanrd <- mean(rd(id))  # mean ret for day id
#   meanrd <- 0            # or meanrd == 0 (no mean correction)
   for (ii in 1:Nrvq)
    {                    # for each rv-interval
      sumr2 <- 0
      for (is in 1:mrv)
       {  tr <- tr + 1;  hhmr <- r[tr] - meanrd; 
          sumr2 <- sumr2 + hhmr*hhmr      # sum of r2 / (r-mu)^2
       }
      trv <- trv + 1
      rv[trv] <- sumr2
    }
   if ( id < 140 )      # for day 140 (no overnight at the end)     
    {
     trv <- trv + 1      # for the overnight vola
     tr  <- tr + 1
     rv[trv] <- (r[tr] - meanrd)*(r[tr] - meanrd)
    }
 }
########
rv <- sqrt(rv)  # realized variance --> realized volatility
                              
dirv_a <- function(d) { (d-1)*Nrvq + (d-1) + 1 }   # begin of day / volatility
dirv_e <- function(d) {  dirv_a(d) + Nrvq -1  }    # end of day  / volatility
dirv_o <- function(d) {  dirv_a(d) + Nrvq }  # overnight vola, end of day, d =/=140
 
rvd <- function(d) { rv[dirv_a(d):dirv_e(d)] }         # without overnight return
#rvd <- function(d) { rv[dirv_a(d):dirv_o(d)] }         # with overnight return
      
rvrv <- as.ts( cbind( rvd(5), rvd(60), rvd(100)) )
plot(rvrv)

###
plot( density(rvd(1)), main="within day rv" )
lines( density(rvd(70)) , col="blue" )
lines( density(rvd(140)) , col="red" )

###########################################################
### Estimation of rv the withinday U-shape by panel methods
# Data layout:
# ... days ... : U
#              : high 
#      1
#      /
#      2
#      h
#      o
#      u       : low
#      r
#
#      u
#      n
#      i
#      t
#      s       : high
#------------- :
# c1 ...   cn  : different constants/level shifts in vola 
##########################################################
# no overnight vola


## EXERCISE:
#  (1) Compare distributions of within day/plus overnight prices.
#  (2) Compare distributions of within day/plus overnight returns.
#  (3) Comapre distributions of within day/plus overnight realized volatility.
#  (3.a) Use different degrees of aggregation for rv.
#  (3.b) Why do we not consider i.g. rv for 1-minute data?