FRE6871_Lecture_7.Rnw

% FRE6871_Lecture_7
% Define knitr options
% !Rnw weave=knitr
% Set global chunk options
<<knitr_setup,include=FALSE,cache=FALSE>>=
library(knitr)
opts_chunk$set(prompt=TRUE, eval=FALSE, tidy=FALSE, strip.white=FALSE, comment=NA, highlight=FALSE, message=FALSE, warning=FALSE, size="tiny", fig.width=4, fig.height=4)
options(width=80, dev="pdf")
options(digits=3)
thm <- knit_theme$get("acid")
knit_theme$set(thm)
@


% Define document options
\documentclass[9pt]{beamer}
\DeclareMathSizes{8pt}{6pt}{6pt}{5pt}
\mode<presentation>
\usetheme{AnnArbor}
% \usecolortheme{whale}
% Uncover everything in a step-wise fashion
% \beamerdefaultoverlayspecification{<+->}
% tikz package for plotting and tables
\usepackage{tikz}
\usetikzlibrary{positioning}
\usepackage{array}
\usepackage{multirow}
% mathtools package for math symbols
% \usepackage{mathtools}
\usepackage{bbm}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage[latin1]{inputenc}
\usepackage{hyperref}
\usepackage{fancybox}
\usepackage{url}
\usepackage[backend=bibtex,style=alphabetic]{biblatex} % bibstyle=numeric
\usepackage{animate}
\usepackage{listings}
\usepackage{xcolor}
\definecolor{anti_flashwhite}{rgb}{0.95, 0.95, 0.96}
\definecolor{cmd_background}{rgb}{0.2, 0.2, 0.0}
\definecolor{vba_background}{rgb}{0.0, 0.0, 0.9}
% \bibliographystyle{amsalpha} % doesn't work
% \addbibresource{FRE_lectures.bib}
% \addbibresource[location=remote]{http://www.citeulike.org/user/jerzyp}
\renewcommand\bibfont{\footnotesize}
\renewcommand{\pgfuseimage}[1]{\scalebox{0.75}{\includegraphics{#1}}} % scale bib icons
\setbeamertemplate{bibliography item}[text] % set bib icons
% \setbeamertemplate{bibliography item}{} % remove bib icons

% \usepackage{enumerate}
% \let\emph\textbf
% \let\alert\textbf
% Define colors for hyperlinks
\definecolor{links}{HTML}{2A1B81}
\hypersetup{colorlinks=true,linkcolor=,urlcolor=links}
% Make url text scriptsize
\renewcommand\UrlFont{\scriptsize}
% Make institute text italic and small
\setbeamerfont{institute}{size=\small,shape=\itshape}
\setbeamerfont{date}{size=\small}
\setbeamerfont{block title}{size=\normalsize} % shape=\itshape
\setbeamerfont{block body}{size=\footnotesize}


% Title page setup
\title[FRE6871 Lecture\#7]{FRE6871 \texttt{R} in Finance}
\subtitle{Lecture\#7, Fall 2024}

\institute[NYU Tandon]{NYU Tandon School of Engineering}
\titlegraphic{\includegraphics[scale=0.2]{image/tandon_long_color.png}}
\author[Jerzy Pawlowski]{Jerzy Pawlowski \emph{\href{mailto:jp3900@nyu.edu}{jp3900@nyu.edu}}}
% \email{jp3900@poly.edu}
\date{October 21, 2024}
% \date{\today}
% \pgfdeclareimage[height=0.5cm]{university-logo}{engineering_long_white}
% \logo{\pgfuseimage{engineering_long_white}}


%%%%%%%%%%%%%%%
\begin{document}


%%%%%%%%%%%%%%%
\maketitle



%%%%%%%%%%%%%%%
\section{Writing and Reading Data from Files}


%%%%%%%%%%%%%%%
\subsection{Writing Text Strings}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{cat()} concatenates strings and writes them to standard output or to files.
      \vskip1ex
      \texttt{cat()} parses its argument character string and its escape sequences (\texttt{"\textbackslash{}"}), but doesn't return a value.
      \vskip1ex
      The function \texttt{print()} doesn't interpret its argument, and simply prints it to standard output and invisibly returns it.
      \vskip1ex
      Typing the name of an object in \texttt{R} implicitly calls \texttt{print()} on that object.
      \vskip1ex
      The function \texttt{save()} writes objects to compressed binary \texttt{.RData} files.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
cat("Enter\ttab")  # Cat() parses backslash escape sequences
print("Enter\ttab")

textv <- print("hello")
textv  # Print() returns its argument

# Create string
textv <- "Title: My Text\nSome numbers: 1,2,3,...\nRprofile files contain code executed at R startup,\n"

cat(textv, file="mytext.txt")  # Write to text file

cat("Title: My Text",  # Write several lines to text file
    "Some numbers: 1,2,3,...",
    "Rprofile files contain code executed at R startup,",
    file="mytext.txt", sep="\n")

save(textv, file="mytext.RData")  # Write to binary file
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Displaying Numeric Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{print()} displays numeric data objects, with the number of digits given by the global option \texttt{"digits"}.
      \vskip1ex
      The function \texttt{sprintf()} returns strings formatted from text strings and numeric data.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=TRUE>>=
print(pi)
print(pi, digits=10)
getOption("digits")
foo <- 12
bar <- "weeks"
sprintf("There are %i %s in the year", foo, bar)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Reading Text from Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{scan()} reads text or data from a file and returns it as a vector or a list.
      \vskip1ex
      The function \texttt{readLines()} reads lines of text from a connection (file or console), and returns them as a vector of \texttt{character} strings.
      \vskip1ex
      The function \texttt{readline()} reads a single line from the console, and returns it as a \texttt{character} string.
      \vskip1ex
      The function \texttt{file.show()} reads text or data from a file and displays in editor.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Read text from file
scan(file="mytext.txt", what=character(), sep="\n")

# Read lines from file
readLines(con="mytext.txt")

# Read text from console
inputv <- readline("Enter a number: ")
class(inputv)
# Coerce to numeric
inputv <- as.numeric(inputv)

# Read text from file and display in editor:
# file.show("mytext.txt")
# file.show("mytext.txt", pager="")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading \protect\emph{Data Frames} from \protect\emph{Text} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The functions \texttt{write.table()} and \texttt{read.table()} write and read \emph{data frames} from text files.
      \vskip1ex
      \texttt{write.table()} coerces objects to \emph{data frames} before it writes them.
      \vskip1ex
      \texttt{read.table()} returns a \emph{data frame}, without coercing non-numeric values to \texttt{factors} (so no need for the option \texttt{stringsAsFactors=FALSE}).
      \vskip1ex
      \texttt{write.table()} and \texttt{read.table()} can be used to write and read matrices from text files, but they have to be coerced back to matrices.
      \vskip1ex
      \texttt{write.table()} and \texttt{read.table()} are inefficient for very large data sets.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
setwd("/Users/jerzy/Develop/lecture_slides/data")
dframe <- data.frame(type=c("rose", "daisy", "tulip"),
  color=c("red", "white", "yellow"),
  price=c(1.5, 0.5, 1.0),
  row.names=c("flower1", "flower2", "flower3"))  # end data.frame
matv <- matrix(sample(1:12), ncol=3,
  dimnames=list(NULL, c("col1", "col2", "col3")))
rownames(matv) <- paste("row", 1:NROW(matv), sep="")
# Write data frame to text file, and then read it back
write.table(dframe, file="florist.txt")
readf <- read.table(file="florist.txt")
readf  # A data frame
all.equal(readf, dframe)

# Write matrix to text file, and then read it back
write.table(matv, file="matrix.txt")
readmat <- read.table(file="matrix.txt")
readmat  # write.table() coerced matrix to data frame
class(readmat)
all.equal(readmat, matv)
# Coerce from data frame back to matrix
readmat <- as.matrix(readmat)
class(readmat)
all.equal(readmat, matv)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Copying \protect\emph{Data Frames} Between the \protect\emph{clipboard} and \texttt{R}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Data frames} stored in the \emph{clipboard} can be copied into \texttt{R} using the function \texttt{read.table()}.
      \vskip1ex
      \emph{Data frames} in \texttt{R} can be copied into the \emph{clipboard} using the function \texttt{write.table()}.
      \vskip1ex
      This allows convenient copying of \emph{data frames} between \texttt{R} and \texttt{Excel}.
      \vskip1ex
      \emph{Data frames} can also be manipulated directly in the \texttt{R} spreadsheet-style data editor.
      \vskip1ex
      Copying and pasting between the \emph{clipboard} and \texttt{R} works well on Windows, but not on MacOS.  There are some workarounds for MacOS:\\
      \href{https://stackoverflow.com/questions/10959521/how-to-write-to-clipboard-on-ubuntu-linux-in-r}{\emph{Copy\_paste\_between\_R\_and\_clipboard}}
    \column{0.5\textwidth}
      \vspace{-1em}
        <<echo=TRUE,eval=FALSE>>=
# Create a data frame
dframe <- data.frame(small=c(3, 5), medium=c(9, 11), large=c(15, 13))

# Launch spreadsheet-style data editor
dframe <- edit(dframe)

# Copy the data frame to clipboard
write.table(x=dframe, file="clipboard", sep="\t")

# Wrapper function for copying data frame from R into clipboard
# by default, data is tab delimited, with a header
write_clip <- function(data, namev=FALSE, col.names=TRUE, ...) {
  write.table(x=data, file="clipboard", sep="\t",
              row.names=namev, col.names=col.names, ...)
}  # end write_clip

write_clip(data=dframe)

# Wrapper function for copying data frame from clipboard into R
# by default, data is tab delimited, with a header
read_clip <- function(file="clipboard", sep="\t", header=TRUE, ...) {
  read.table(file=file, sep=sep, header=header, ...)
}  # end read_clip

dframe <- read.table("clipboard", header=TRUE)
dframe <- read_clip()
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading \protect\emph{Data Frames} From \texttt{.csv} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The easiest way to share data between \texttt{R} and \texttt{Excel} is through \texttt{.csv} files.
      \vskip1ex
      The functions \texttt{write.csv()} and \texttt{read.csv()} write and read \emph{data frames} from \texttt{.csv} format files.
      \vskip1ex
      The functions \texttt{write.csv()} and \texttt{read.csv()} write and read \emph{data frames} from \texttt{.csv} format files.
      \vskip1ex
      These functions are \emph{wrappers} for \texttt{write.table()} and \texttt{read.table()}.
      \vskip1ex
      \texttt{read.csv()} doesn't coerce non-numeric values to \texttt{factors}, so no need for the option \texttt{stringsAsFactors=FALSE}.
      \vskip1ex
      \texttt{read.csv()} reads row names as an extra column, unless the \texttt{row.names=1} argument is used.
      \vskip1ex
      The argument \texttt{"row.names"} accepts either the number or the name of the column containing the row names.
      \vskip1ex
      The \texttt{*.csv()} functions are very inefficient for large data sets.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Write data frame to CSV file, and then read it back
write.csv(dframe, file="florist.csv")
readf <- read.csv(file="florist.csv")
readf  # the row names are read in as extra column
# Restore row names
rownames(readf) <- readf[, 1]
readf <- readf[, -1]  # Remove extra column
readf
all.equal(readf, dframe)
# Read data frame, with row names from first column
readf <- read.csv(file="florist.csv", row.names=1)
readf
all.equal(readf, dframe)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading \protect\emph{Data Frames} From \texttt{.csv} Files (cont.)}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The functions \texttt{write.csv()} and \texttt{read.csv()} can write and read \emph{data frames} from \texttt{.csv} format files \emph{without using row names}.
      \vskip1ex
      Row names can be omitted from the output file by calling \texttt{write.csv()} with the argument \texttt{row.names=FALSE}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Write data frame to CSV file, without row names
write.csv(dframe, row.names=FALSE, file="florist.csv")
readf <- read.csv(file="florist.csv")
readf  # A data frame without row names
all.equal(readf, dframe)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Reading Data From Very Large \texttt{.csv} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Data from very large \texttt{.csv} files can be read in small chunks instead of all at once.
      \vskip1ex
      The function \texttt{file()} opens a connection to a file or an internet website \texttt{URL}.
      \vskip1ex
      The function \texttt{read.csv()} with the argument \texttt{"nrows"} reads only the specified number of rows from a connection and returns a \emph{data frame}.  The connection pointer is reset to the next row.
      \vskip1ex
      The function \texttt{read.csv()} with the argument \texttt{"nrows"} allows reading data sequentially from very large files that wouldn't fit into memory.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Open a read connection to a file
filecon = file("/Users/jerzy/Develop/lecture_slides/data/etf_prices_crsp.csv", "r")
# Read the first 10 rows
data10 <- read.csv(filecon, nrows=10)
# Read another 10 rows
data20 <- read.csv(filecon, nrows=10, header=FALSE)
colnames(data20) <- colnames(data10)
# Close the connection to the file
close(filecon)
# Open a read connection to a file
filecon = file("/Users/jerzy/Develop/lecture_slides/data/etf_prices_crsp.csv", "r")
# Read the first 1000 rows
data10 <- read.csv(filecon, nrows=1e3)
colv <- colnames(data10)
# Write to a file
countv <- 1
write.csv(data10, paste0("/Users/jerzy/Develop/data/temp/etf_prices_", countv, ".csv"))
# Read remaining rows in a loop 10 rows at a time
# Can produce error without getting to end of file
while (isOpen(filecon)) {
  datav <- read.csv(filecon, nrows=1e3)
  colnames(datav) <- colv
  write.csv(datav, paste0("/Users/jerzy/Develop/data/temp/etf_prices_", countv, ".csv"))
  countv <- countv + 1
}  # end while
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Matrices From \texttt{.csv} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The functions \texttt{write.csv()} and \texttt{read.csv()} can write and read matrices from \texttt{.csv} format files.
      \vskip1ex
      If row names can be omitted in the output file, then \texttt{write.csv()} can be called with argument \texttt{row.names=FALSE}.
      \vskip1ex
      If the input file doesn't contain row names, then \texttt{read.csv()} can be called without the \texttt{"row.names"} argument.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Write matrix to csv file, and then read it back
write.csv(matv, file="matrix.csv")
readmat <- read.csv(file="matrix.csv", row.names=1)
readmat  # Read.csv() reads matrix as data frame
class(readmat)
readmat <- as.matrix(readmat)  # Coerce to matrix
all.equal(readmat, matv)
write.csv(matv, row.names=FALSE,
          file="matrix_ex_rows.csv")
readmat <- read.csv(file="matrix_ex_rows.csv")
readmat <- as.matrix(readmat)
readmat  # A matrix without row names
all.equal(readmat, matv)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Matrices (cont.)}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      There are several ways of writing and reading matrices from \texttt{.csv} files, with tradeoffs between simplicity, data size, and speed.
      \vskip1ex
      The function \texttt{write.matrix()} writes a matrix to a text file, without its row names.
      \vskip1ex
      \texttt{write.matrix()} is part of package \emph{MASS}.
      \vskip1ex
      The advantage of function \texttt{scan()} is its speed, but it doesn't handle row names easily.
      \vskip1ex
      Removing row names simplifies the writing and reading of matrices.
      \vskip1ex
      The function \texttt{readLines} reads whole lines and returns them as single strings.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
setwd("/Users/jerzy/Develop/lecture_slides/data")
library(MASS)  # Load package "MASS"
# Write to CSV file by row - it's very SLOW!!!
MASS::write.matrix(matv, file="matrix.csv", sep=",")
# Read using scan() and skip first line with colnames
readmat <- scan(file="matrix.csv", sep=",", skip=1,
  what=numeric())
# Read colnames
colv <- readLines(con="matrix.csv", n=1)
colv  # this is a string!
# Convert to char vector
colv <- strsplit(colv, split=",")[[1]]
readmat  # readmat is a vector, not matrix!
# Coerce by row to matrix
readmat <- matrix(readmat, ncol=NROW(colv), byrow=TRUE)
# Restore colnames
colnames(readmat) <- colv
readmat
# Scan() is a little faster than read.csv()
library(microbenchmark)
summary(microbenchmark(
  read_csv=read.csv("matrix.csv"),
  scan=scan(file="matrix.csv", sep=",",
    skip=1, what=numeric()),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Reading Matrices Containing Bad Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Very often data that is read from external sources contains elements with bad data.
      \vskip1ex
      An example of bad data are \texttt{character} strings within sets of \texttt{numeric} data.
      \vskip1ex
      Columns of numeric data that contain strings are coerced to \texttt{character} or \texttt{factor}, when they're read by \texttt{read.csv()}.
      \vskip1ex
      The function \texttt{as.numeric()} coerces complex data objects into \texttt{numeric} vectors, and removes all their \emph{attributes}.
      \vskip1ex
      \texttt{as.numeric()} coerces strings that don't represent numbers into \texttt{NA} values.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Read data from a csv file, including row names
matv <- read.csv(file="data/matrix_bad.csv", row.names=1)
matv
class(matv)
# Columns with bad data are character or factor
sapply(matv, class)
# Coerce character column to numeric
matv$col2 <- as.numeric(matv$col2)
# Or
# Copy row names
namev <- rownames(matv)
# sapply loop over columns and coerce to numeric
matv <- sapply(matv, as.numeric)
# Restore row names
rownames(matv) <- namev
# Replace NAs with zero
matv[is.na(matv)] <- 0
# matrix without NAs
matv
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Time Series From \protect\emph{Text} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package \emph{zoo} contains functions \texttt{write.zoo()} and \texttt{read.zoo()} for writing and reading \emph{zoo} time series from \texttt{.txt} and \texttt{.csv} files.
      \vskip1ex
      The functions \texttt{write.zoo()} and \texttt{read.zoo()} are \emph{wrappers} for \texttt{write.table()} and \texttt{read.table()}.
      \vskip1ex
      The function \texttt{write.zoo()} writes the \emph{zoo} series index as a character string in quotations \texttt{""}, to make it easier to read (parse) by \texttt{read.zoo()}.
      \vskip1ex
      Users may also directly use \texttt{write.table()} and \texttt{read.table()}, instead of \texttt{write.zoo()} and \texttt{read.zoo()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:4)),eval=FALSE>>=
setwd("/Users/jerzy/Develop/lecture_slides/data")
rm(list=ls())
# Initialize the random number generator
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
library(zoo)  # Load package zoo
# Create zoo with Date index
datev <- seq(from=as.Date("2013-06-15"), by="day",
              length.out=100)
pricev <- zoo(rnorm(NROW(datev)), order.by=datev)
head(pricev, 3)
# Write zoo series to text file, and then read it back
write.zoo(pricev, file="pricev.txt")
pricezoo <- read.zoo("pricev.txt")  # Read it back
all.equal(pricezoo, pricev)
# Perform the same using write.table() and read.table()
# First coerce pricev into data frame
dframe <- as.data.frame(pricev)
dframe <- cbind(datev, dframe)
# Write pricev to text file using write.table
write.table(dframe, file="pricev.txt",
            row.names=FALSE, col.names=FALSE)
# Read data frame from file
pricezoo <- read.table(file="pricev.txt")
sapply(pricezoo, class)  # A data frame
# Coerce data frame into pricev
pricezoo <- zoo::zoo(
  drop(as.matrix(pricezoo[, -1])),
  order.by=as.Date(pricezoo[, 1]))
all.equal(pricezoo, pricev)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Time Series From \texttt{.csv} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      By default the functions \texttt{zoo::write.zoo()} and \texttt{zoo::read.zoo()} write data in \emph{space}-delimited text format, but they can also write to \emph{comma}-delimited \texttt{.csv} files by passing the parameter \texttt{sep=","}.
      \vskip1ex
      Single column \emph{zoo} time series usually don't have a dimension attribute, and they don't have a column name, unlike multi-column \emph{zoo} time series, and this can cause hard to detect bugs.
      \vskip1ex
      It's best to always pass the argument \texttt{"col.names=TRUE"} to the function \texttt{write.zoo()}, to make sure it writes a column name for a single column \emph{zoo} time series.
      \vskip1ex
      Reading a \texttt{.csv} file containing a single column of data using the function \texttt{read.zoo()} produces a \emph{zoo} time series with a \texttt{NULL} dimension, unless the argument \texttt{"drop=FALSE"} is passed to \texttt{read.zoo()}.
      \vskip1ex
      Users may also directly use \texttt{write.table()} and \texttt{read.table()}, instead of \texttt{write.zoo()} and \texttt{read.zoo()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(zoo)  # Load package zoo
# Write zoo series to CSV file, and then read it back
write.zoo(pricev, file="zooseries.csv", sep=",", col.names=TRUE)
pricezoo <- read.zoo(file="zooseries.csv",
  header=TRUE, sep=",", drop=FALSE)
all.equal(pricev, drop(pricezoo))
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Time Series With \protect\emph{Date-time} Index}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{write.zoo()} writes \emph{zoo} time series into \texttt{.csv} files, but it doesn't format the time at midnight properly.
      \vskip1ex
      The function \texttt{write.table()} writes \emph{zoo} time series into \texttt{.csv} files, and it formats the time at midnight properly.
      \vskip1ex
      If the index of a \emph{zoo} time series is a \emph{date-time}, then \texttt{write.zoo()} writes the date and time fields as character strings separated by a \emph{space} between them, inside quotations \texttt{""}.
      \vskip1ex
      The functions \texttt{read.csv.zoo()} and \texttt{read.zoo()} read \emph{zoo} time series from \texttt{.csv} files.
      \vskip1ex
      Very often \texttt{.csv} files contain custom \emph{date-time} formats, which need to be passed as parameters into \texttt{read.zoo()} for proper formatting.
      \vskip1ex
      The \texttt{"FUN"} argument of \texttt{read.zoo()} accepts a function for coercing the date and time columns of the input data into a \emph{date-time} object suitable for the \emph{zoo} index.
      \vskip1ex
      The function \texttt{as.POSIXct()} coerces \texttt{character} strings into \texttt{POSIXct} \emph{date-time} objects.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
# Initialize the random number generator
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
# Create zoo with POSIXct date-time index
datev <- seq(from=as.POSIXct("2014-07-14"),
              by="hour", length.out=100)
zooseries <- zoo(rnorm(NROW(datev)), order.by=datev)
head(zooseries, 3)
# Write zoo series to CSV file using write.zoo()
write.zoo(zooseries, file="zooseries.csv", sep=",", col.names=TRUE)
# Read from CSV file using read.csv.zoo() - doesn't work
zooread <- read.csv.zoo(file="zooseries.csv", header=FALSE,
  format="%Y-%m-%d %H:%M:%S", tz="America/New_York")
# Read from CSV file using read.zoo() - error
zooread <- read.zoo(file="zooseries.csv", header=FALSE, 
  sep=",", FUN=as.POSIXct, format="%Y-%m-%d %H:%M:%S")
# Write zoo series to CSV file using write.table()
write.table(zooseries, file="zooseries.csv", sep=",", 
            row.names=TRUE, col.names=FALSE)
# Read from CSV file using read.zoo() with format argument
zooread <- read.zoo(file="zooseries.csv", header=FALSE, 
  sep=",", FUN=as.POSIXct, format="%Y-%m-%d %H:%M:%S")
all.equal(zooseries, zooread) # Works
# Coerce zoo series into data frame with custom date format
dframe <- as.data.frame(zooseries)
rownames(dframe) <- format(index(zooseries), format="%m-%d-%Y %H:%M:%S")
# Write zoo series to csv file using write.table
write.table(dframe, file="zooseries.csv", sep=",", 
            row.names=TRUE, col.names=FALSE)
# Read from CSV file using read.zoo()
zooread <- read.zoo(file="zooseries.csv", header=FALSE, sep=",", 
  FUN=as.POSIXct, format="%m-%d-%Y %H:%M:%S")
all.equal(zooseries, zooread) # Works
# Or using read.csv.zoo()
zooread <- read.csv.zoo(file="zooseries.csv", header=FALSE,
  format="%m-%d-%Y %H:%M:%S", tz="America/New_York")
head(zooread, 3)
all.equal(zooseries, zooread, check.attributes=FALSE) # Works
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Reading Time Series With \texttt{Numeric} \protect\emph{Date-time} Index}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      If the index of a time series is \texttt{numeric} (representing the \emph{moment of time}, either as the number of days or seconds), then it must be coerced to a proper \emph{date-time} class.
      \vskip1ex
      A convenient way of reading time series with a numeric index is by using \texttt{read.table()}, and then coercing the \emph{data frame} into a time series.
      \vskip1ex
      The function \texttt{as.POSIXct.numeric()} coerces a \texttt{numeric} value representing the \emph{moment of time} into a \texttt{POSIXct} \emph{date-time}, equal to the \emph{clock time} in the local \emph{time zone}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Read time series from CSV file, with numeric date-time
datazoo <- read.table(file="/Users/jerzy/Develop/lecture_slides/data/es_ohlc.csv",
  header=TRUE, sep=",")
# A data frame
class(datazoo)
sapply(datazoo, class)
# Coerce data frame into xts series
datazoo <- xts::xts(as.matrix(datazoo[, -1]),
  order.by=as.POSIXct.numeric(datazoo[, 1], tz="America/New_York",
                              origin="1970-01-01"))
# An xts series
class(datazoo)
head(datazoo, 3)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Passing Arguments to the \texttt{save()} Function}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{save()} writes objects to a binary file.
      \vskip1ex
      Object names can be passed into \texttt{save()} either through the \texttt{"..."} argument, or the \texttt{"list"} argument.
      \vskip1ex
      Objects passed through the \texttt{"..."} argument are not evaluated, so they must be either object names or character strings.
      \vskip1ex
      Object names aren't surrounded by quotes \texttt{""}, while character strings that represent object names are surrounded by quotes \texttt{""}.
      \vskip1ex
      Objects passed through the \texttt{"list"} argument are evaluated, so they may be variables containing character strings.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-1)>>=
rm(list=ls())  # Delete all objects in workspace
var1 <- 1; var2 <- 2
ls()  # List all objects
ls()[1]  # List first object
args(save)  # List arguments of save function
# Save "var1" to a binary file using string argument
save("var1", file="my_data.RData")
# Save "var1" to a binary file using object name
save(var1, file="my_data.RData")
# Save multiple objects
save(var1, var2, file="my_data.RData")
# Save first object in list by passing to "..." argument
# ls()[1] is not evaluated
save(ls()[1], file="my_data.RData")
# Save first object in list by passing to "list" argument
save(list=ls()[1], file="my_data.RData")
# Save whole list by passing it to the "list" argument
save(list=ls(), file="my_data.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Lists of Objects}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{load()} reads data from \texttt{.RData} files, and \emph{invisibly} returns a vector of names of objects created in the workspace.
      \vskip1ex
      The vector of names can be used to manipulate the objects in loops, or to pass them to functions.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
rm(list=ls())  # Delete all objects in workspace
# Load objects from file
loadobj <- load(file="my_data.RData")
loadobj  # vector of loaded objects
ls()  # List objects
# Assign new values to objects in  global environment
sapply(loadobj, function(symboln) {
  assign(symboln, runif(1), envir=globalenv())
})  # end sapply
ls()  # List objects
# Assign new values to objects using for loop
for (symboln in loadobj) {
  assign(symboln, runif(1))
}  # end for
ls()  # List objects
# Save vector of objects
save(list=loadobj, file="my_data.RData")
# Remove only loaded objects
rm(list=loadobj)
# Remove the object "loadobj"
rm(loadobj)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Saving Output of \texttt{R} to a File}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{sink()} diverts \texttt{R} \emph{text} output (excluding graphics) to a file, or ends the diversion.
      \vskip1ex
      Remember to call \texttt{sink()} to end the diversion!
      \vskip1ex
      The function \texttt{pdf()} diverts graphics output to a \emph{pdf} file (text output isn't diverted), in vector graphics format.
      \vskip1ex
      The functions \texttt{png()}, \texttt{jpeg()}, \texttt{bmp()}, and \texttt{tiff()} divert graphics output to graphics files (text output isn't diverted).
      \vskip1ex
      The function \texttt{dev.off()} ends the diversion.
    \column{0.5\textwidth}
      \vspace{-1em}
        <<echo=TRUE,eval=FALSE>>=
sink("sinkdata.txt")# Redirect text output to file

cat("Redirect text output from R\n")
print(runif(10))
cat("\nEnd data\nbye\n")

sink()  # turn redirect off

pdf("Rgraph.pdf", width=7, height=4)  # Redirect graphics to pdf file

cat("Redirect data from R into pdf file\n")
myvar <- seq(-2*pi, 2*pi, len=100)
plot(x=myvar, y=sin(myvar), main="Sine wave",
   xlab="", ylab="", type="l", lwd=2, col="red")
cat("\nEnd data\nbye\n")

dev.off()  # turn pdf output off

png("r_plot.png")  # Redirect graphics output to png file

cat("Redirect graphics from R into png file\n")
plot(x=myvar, y=sin(myvar), main="Sine wave",
 xlab="", ylab="", type="l", lwd=2, col="red")
cat("\nEnd data\nbye\n")

dev.off()  # turn png output off
      @
  \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{High Performance Data Management}


%%%%%%%%%%%%%%%
\subsection{Package \protect\emph{data.table} for High Performance Data Management}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package
      \href{https://cran.r-project.org/web/packages/data.table/}{\emph{data.table}}
      is designed for high performance data management.
      \vskip1ex
      The package \emph{data.table} implements \emph{data table} objects, which are a special type of \emph{data frame}, and an extension of the \emph{data frame} class.
      \vskip1ex
      \emph{Data tables} are faster and more convenient to work with than \emph{data frames}.
      \vskip1ex
      \emph{data.table} functions are optimized for high performance (speed), because they are written in \texttt{C++} and they perform operations by reference (in place), without copying data in memory.
      \vskip1ex
      Some of the attractive features of package \emph{data.table} are:
      \begin{itemize}
        \item Syntax is analogous to SQL,
        \item Very fast writing and reading from files,
        \item Very fast sorting and merging operations,
        \item Subsetting using multiple logical clauses,
        \item Columns of type \texttt{character} are never converted to factors,
      \end{itemize}
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Install package data.table
install.packages("data.table")
# Load package data.table
library(data.table)
# Get documentation for package data.table
# Get short description
packageDescription("data.table")
# Load help page
help(package="data.table")
# List all datasets in "data.table"
data(package="data.table")
# List all objects in "data.table"
ls("package:data.table")
# Remove data.table from search path
detach("package:data.table")
      @
      The package \emph{data.table} has extensive documentation:\\
      \hskip1em\url{https://cran.r-project.org/web/packages/data.table/vignettes/datatable-intro.html}\\
      \hskip1em\url{https://github.com/Rdatatable/data.table/wiki}
      \vskip1ex
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{\protect\emph{Data Table} Objects}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Data table} objects are a special type of \emph{data frame}, and are derived from the class \texttt{data.frame}.
      \vskip1ex
      \emph{Data table} objects resemble databases, with columns of different types of data, and rows of records containing individual observations.
      \vskip1ex
      The function \texttt{data.table::data.table()} creates a \emph{data table} object.
      \vskip1ex
      \emph{Data table} columns can be referenced directly by their names (without quotes), and their rows can be referenced without a following comma.
      \vskip1ex
      When a \emph{data table} is printed (by typing its name) then only the top \texttt{5} and bottom \texttt{5} rows are displayed (unless \texttt{getOption("datatable.print.nrows")} is less than \texttt{100}).
      \vskip1ex
      The operator \texttt{.N} returns the number of observations (rows) in the \emph{data table}.
      \vskip1ex
      \emph{Data table} computations are usually much faster than equivalent \texttt{R} computations, but not always.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Create a data table
library(data.table)
dtable <- data.table::data.table(
  col1=sample(7), col2=sample(7), col3=sample(7))
# Print dtable
class(dtable); dtable
# Column referenced without quotes
dtable[, col2]
# Row referenced without a following comma
dtable[2]
# Print option "datatable.print.nrows"
getOption("datatable.print.nrows")
options(datatable.print.nrows=10)
getOption("datatable.print.nrows")
# Number of rows in dtable
NROW(dtable)
# Or
dtable[, NROW(col1)]
# Or
dtable[, .N]
# microbenchmark speed of data.table syntax
library(microbenchmark)
summary(microbenchmark(
  dt=dtable[, .N],
  rcode=NROW(dtable),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Data Using Package \protect\emph{data.table}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The easiest way to share data between \texttt{R} and \texttt{Excel} is through \texttt{.csv} files.
      \vskip1ex
      The function \texttt{data.table::fread()} reads from \texttt{.csv} files and returns a \emph{data table} object of class \texttt{data.table}.
      \vskip1ex
      \emph{Data table} objects are a special type of \emph{data frame}, and are derived from the class \texttt{data.frame}.
      \vskip1ex
      The function \texttt{data.table::fread()} is over \texttt{6} times faster than \texttt{read.csv()}!
      \vskip1ex
      The function \texttt{data.table::fwrite()} writes to \texttt{.csv} files over \texttt{12} times faster than the function \texttt{write.csv()}, and \texttt{300} times faster than function \texttt{cat()}!
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Read a data table from CSV file
dirn <- "/Users/jerzy/Develop/lecture_slides/data/"
filen <- file.path(dirn, "weather_delays14.csv")
dtable <- data.table::fread(filen)
class(dtable); dim(dtable)
dtable
# fread() reads the same data as read.csv()
all.equal(read.csv(filen),
          setDF(data.table::fread(filen)))
# fread() is much faster than read.csv()
library(microbenchmark)
summary(microbenchmark(
  rcode=read.csv(filen),
  fread=setDF(data.table::fread(filen)),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
# Write data table to file in different ways
data.table::fwrite(dtable, file="dtable.csv")
write.csv(dtable, file="dtable2.csv")
cat(unlist(dtable), file="dtable3.csv")
# microbenchmark speed of data.table::fwrite()
summary(microbenchmark(
  fwrite=data.table::fwrite(dtable, file="dtable.csv"),
  write_csv=write.csv(dtable, file="dtable2.csv"),
  cat=cat(unlist(dtable), file="dtable3.csv"),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Subsetting \protect\emph{Data Table} Objects}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The square braces (brackets) \texttt{"[]"} operator subsets (references) the rows and columns of \emph{data tables}.
      \vskip1ex
      \emph{Data table} rows can be subset without a following comma.
      \vskip1ex
      \emph{Data table} columns can be referenced directly by their names (without quotes, as if they were variables), after a comma.
      \vskip1ex
      Multiple \emph{data table} columns can be referenced by passing a list of names.
      \vskip1ex
      The brackets \texttt{"[]"} operator is a \emph{data.table} function, and all the commands inside the brackets \texttt{"[]"} are executed using code from the package \emph{data.table}.
      \vskip1ex
      The dot \texttt{.()} operator is equivalent to the list function \texttt{list()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Select first five rows of dtable
dtable[1:5]
# Select rows with JFK flights
jfkf <- dtable[origin=="JFK"]
# Select rows JFK flights in June
jfkf <- dtable[origin=="JFK" & month==6]
# Select rows without JFK flights
jfkf <- dtable[!(origin=="JFK")]
# Select flights with carrier_delay
dtable[carrier_delay > 0]
# Select column of dtable and return a vector
head(dtable[, origin])
# Select column of dtable and return a dtable, not vector
head(dtable[, list(origin)])
head(dtable[, .(origin)])
# Select two columns of dtable
dtable[, list(origin, month)]
dtable[, .(origin, month)]
columnv <- c("origin", "month")
dtable[, ..columnv]
dtable[, month, origin]
# Select two columns and rename them
dtable[, .(orig=origin, mon=month)]
# Select all columns except origin
head(dtable[, !"origin"])
head(dtable[, -"origin"])
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Performing Computations on \protect\emph{Data Table} Columns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      If the second argument in the brackets \texttt{"[]"} operator is a function of the columns, then the brackets return the result of the function's computations on those columns.
      \vskip1ex
      The second argument in the brackets \texttt{"[]"} can also be a list of functions, in which case the brackets return a vector of computations.
      \vskip1ex
      The brackets \texttt{"[]"} can evaluate most standard \texttt{R} functions, but they are executed using \emph{data.table} code, which is usually much faster than the equivalent \texttt{R} functions.
      \vskip1ex
      The operator \texttt{.N} returns the number of observations (rows) in the \emph{data table}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Select flights with positive carrier_delay
dtable[carrier_delay > 0]
# Number of flights with carrier_delay
dtable[, sum(carrier_delay > 0)]
# Or standard R commands
sum(dtable[, carrier_delay > 0])
# microbenchmark speed of data.table syntax
summary(microbenchmark(
  dt=dtable[, sum(carrier_delay > 0)],
  rcode=sum(dtable[, carrier_delay > 0]),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
# Average carrier_delay
dtable[, mean(carrier_delay)]
# Average carrier_delay and aircraft_delay
dtable[, .(carrier=mean(carrier_delay),
           aircraft=mean(aircraft_delay))]
# Average aircraft_delay from JFK
dtable[origin=="JFK", mean(aircraft_delay)]
# Number of flights from JFK
dtable[origin=="JFK", NROW(aircraft_delay)]
# Or
dtable[origin=="JFK", .N]
# In R
sum(dtable[, origin]=="JFK")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Grouping \protect\emph{Data Table} Computations by Factor Columns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{data table} brackets \texttt{"[]"} operator can accept three arguments: \texttt{[i, j, by]}
      \begin{itemize}
        \item \texttt{i}: the row index to select,
        \item \texttt{j}: a list of columns or functions on columns,
        \item \texttt{by}: the columns of factors to aggregate over.
      \end{itemize}
      The \emph{data table} columns can be \emph{aggregated} over categories (factors) defined by one or more columns passed to the \texttt{"by"} argument.
      \vskip1ex
      The \texttt{"keyby"} argument is similar to \texttt{"by"}, but it sorts the output according to the categories used to group by.
      \vskip1ex
      Multiple \emph{data table} columns can be referenced by passing a list of names.
      \vskip1ex
      The dot \texttt{.()} operator is equivalent to the list function \texttt{list()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Number of flights from each airport
dtable[, .N, by=origin]
# Same, but add names to output
dtable[, .(flights=.N), by=.(airport=origin)]
# Number of AA flights from each airport
dtable[carrier=="AA", .(flights=.N), by=.(airport=origin)]
# Number of flights from each airport and airline
dtable[, .(flights=.N), by=.(airport=origin, airline=carrier)]
# Average aircraft_delay
dtable[, mean(aircraft_delay)]
# Average aircraft_delay from JFK
dtable[origin=="JFK", mean(aircraft_delay)]
# Average aircraft_delay from each airport
dtable[, .(delay=mean(aircraft_delay)), by=.(airport=origin)]
# Average and max delays from each airport and month
dtable[, .(mean_delay=mean(aircraft_delay), max_delay=max(aircraft_delay)),
           by=.(airport=origin, month=month)]
# Average and max delays from each airport and month
dtable[, .(mean_delay=mean(aircraft_delay), max_delay=max(aircraft_delay)),
           keyby=.(airport=origin, month=month)]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Sorting \protect\emph{Data Table} Rows by Columns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Standard \texttt{R} functions can be used inside the brackets \texttt{"[]"} operator.
      \vskip1ex
      The function \texttt{order()} calculates the permutation index, to sort a given vector into ascending order.
      \vskip1ex
      The function \texttt{setorder()} sorts the rows of a \emph{data table} by reference (in place), without copying data in memory.
      \vskip1ex
      \texttt{setorder()} is over \texttt{10} times faster than \texttt{order()}, because it doesn't copy data in memory.
      \vskip1ex
      Several brackets \texttt{"[]"} operators can be chained together to perform several consecutive computations.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Sort ascending by origin, then descending by dest
dtables <- dtable[order(origin, -dest)]
dtables
# Doesn't work outside dtable
order(origin, -dest)
# Sort dtable by reference
setorder(dtable, origin, -dest)
all.equal(dtable, dtables)
# setorder() is much faster than order()
summary(microbenchmark(
  order=dtable[order(origin, -dest)],
  setorder=setorder(dtable, origin, -dest),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
# Average aircraft_delay by month
dtables[, .(mean_delay=mean(aircraft_delay)),
            by=.(month=month)]
# Chained brackets to sort output by month
dtables[, .(mean_delay=mean(aircraft_delay)),
        by=.(month=month)][order(month)]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Subsetting, Computing, and Grouping \protect\emph{Data Table} Objects}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The special symbol \texttt{.SD} selects a subset of a \emph{data table}.
      \vskip1ex
      The symbol \texttt{.SDcols} specifies the columns to select by the symbol \texttt{.SD}.
      \vskip1ex
      Inside the brackets \texttt{"[]"} operator, the \texttt{.SD} symbol can be treated as a virtual \emph{data table}, and standard \texttt{R} functions can be applied to it.
      \vskip1ex
      The \texttt{"by"} argument can be used to group the outputs produced by the functions applied to the \texttt{.SD} symbol.
      \vskip1ex
      If the symbol \texttt{.SDcols} is not defined, then the symbol \texttt{.SD} returns the remaining columns not passed to the \texttt{"by"} operator.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Select weather_delay and aircraft_delay in two different ways
dtable[1:7, .SD,
           .SDcols=c("weather_delay", "aircraft_delay")]
dtable[1:7, .(weather_delay, aircraft_delay)]
# Calculate mean of weather_delay and aircraft_delay
dtable[, sapply(.SD, mean),
           .SDcols=c("weather_delay", "aircraft_delay")]
sapply(dtable[, .SD,
           .SDcols=c("weather_delay", "aircraft_delay")], mean)
# Return origin and dest, then all other columns
dtable[1:7, .SD, by=.(origin, dest)]
# Return origin and dest, then weather_delay and aircraft_delay columns
dtable[1:7, .SD, by=.(origin, dest),
           .SDcols=c("weather_delay", "aircraft_delay")]
# Return first two rows from each month
dtable[, head(.SD, 2), by=.(month)]
dtable[, head(.SD, 2), by=.(month),
           .SDcols=c("weather_delay", "aircraft_delay")]
# Calculate mean of weather_delay and aircraft_delay, grouped by origin
dtable[, lapply(.SD, mean),
           by=.(origin),
           .SDcols=c("weather_delay", "aircraft_delay")]
# Or simply
dtable[, .(weather_delay=mean(weather_delay),
               aircraft_delay=mean(aircraft_delay)),
           by=.(origin)]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Modifying \protect\emph{Data Table} Objects by Reference}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The special assignment operator \texttt{":="} allows modifying \emph{data table} columns by reference (in place), without copying data in memory.
      \vskip1ex
      The computations on columns by reference can be \emph{grouped} over categories defined by one or more columns passed to the \texttt{"by"} argument.
      \vskip1ex
      The computations are recycled to fit the size of each group.
      \vskip1ex
      The selected parts of columns can also be modified by reference, by combining the \texttt{i} and \texttt{j} arguments.
      \vskip1ex
      The special symbols \texttt{.SD} and \texttt{.SDcols} can be used to perform computations on several columns.
      \vskip1ex
      Modifying by reference is several times faster than standard \texttt{R} assignment.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Add tot_delay column
dtable[, tot_delay := (carrier_delay + aircraft_delay)]
head(dtable, 4)
# Delete tot_delay column
dtable[, tot_delay := NULL]
# Add max_delay column grouped by origin and dest
dtable[, max_delay := max(aircraft_delay), by=.(origin, dest)]
dtable[, max_delay := NULL]
# Add date and tot_delay columns
dtable[, c("date", "tot_delay") :=
             list(paste(month, day, year, sep="/"),
                  (carrier_delay + aircraft_delay))]
# Modify select rows of tot_delay column
dtable[month == 12, tot_delay := carrier_delay]
dtable[, c("date", "tot_delay") := NULL]
# Add several columns
dtable[, c("max_carrier", "max_aircraft") := lapply(.SD, max),
       by=.(origin, dest),
       .SDcols=c("carrier_delay", "aircraft_delay")]
# Remove columns
dtable[, c("max_carrier", "max_aircraft") := NULL]
# Modifying by reference is much faster than standard R
summary(microbenchmark(
  dt=dtable[, tot_delay := (carrier_delay + aircraft_delay)],
  rcode=(dtable[, "tot_delay"] <- dtable[, "carrier_delay"] + dtable[, "aircraft_delay"]),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Adding \protect\emph{keys} to \protect\emph{Data Tables} for Fast Binary Search}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{key} of a \emph{data table} is analogous to the row indices of a \emph{data frame}, and it determines the ordering of its rows.
      \vskip1ex
      The function \texttt{data.table::setkey()} adds a \emph{key} to a \emph{data table}, and sorts the \emph{data table} rows by reference according to the key.
      \vskip1ex
      \texttt{setkey()} creates the \emph{key} from one or more columns of the \emph{data frame}.
      \vskip1ex
      Subsetting rows using a \emph{key} can be several times faster than standard \texttt{R}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Add a key based on the "origin" column
setkey(dtable, origin)
haskey(dtable)
key(dtable)
# Select rows with LGA using the key
dtable["LGA"]
all.equal(dtable["LGA"], dtable[origin == "LGA"])
# Select rows with LGA and JFK using the key
dtable[c("LGA", "JFK")]
# Add a key based on the "origin" and "dest" columns
setkey(dtable, origin, dest)
key(dtable)
# Select rows with origin from JFK and MIA
dtable[c("JFK", "MIA")]
# Select rows with origin from JFK and dest to MIA
dtable[.("JFK", "MIA")]
all.equal(dtable[.("JFK", "MIA")],
          dtable[origin == "JFK" & dest == "MIA"])
# Selecting rows using a key is much faster than standard R
summary(microbenchmark(
  with_key=dtable[.("JFK", "MIA")],
  standard_r=dtable[origin == "JFK" & dest == "MIA"],
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Coercing \protect\emph{Data Table} Objects Into \protect\emph{Data Frames}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The functions \texttt{data.table::setDT()} and \texttt{data.table::setDF()} coerce \emph{data frames} to \emph{data tables}, and vice versa.
      \vskip1ex
      The \emph{set} functions \texttt{data.table::set*()} perform their operations by reference (in place), without returning any values or copying data to a new memory location, which makes them very fast.
      \vskip1ex
      \emph{Data table} objects can also be coerced into \emph{data frames} using the function \texttt{as.data.frame()}, but it's much slower because it makes copies of data.
      \vskip1ex
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Create data frame and coerce it to data table
dtable <- data.frame(col1=sample(7), col2=sample(7), col3=sample(7))
class(dtable); dtable
data.table::setDT(dtable)
class(dtable); dtable
# Coerce dtable into data frame
data.table::setDF(dtable)
class(dtable); dtable
# Or
dtable <- data.table:::as.data.frame.data.table(dtable)
# SetDF() is much faster than as.data.frame()
summary(microbenchmark(
  asdataframe=data.table:::as.data.frame.data.table(dtable),
  setDF=data.table::setDF(dtable),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Coercing \texttt{xts} Time Series Into \protect\emph{Data Tables}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      An \texttt{xts} time series can be coerced into a \emph{data table} by first coercing it into a \emph{data frame} and then into a \emph{data table} using the function \texttt{data.table::setDT()}.
      \vskip1ex
      But then the time index of the \texttt{xts} series is coerced into strings, not dates.
      \vskip1ex
      An \texttt{xts} time series can also be coerced directly into a \emph{data table} using the function \texttt{data.table::as.data.table()}.
      \vskip1ex
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Coerce xts to a data frame
pricev <- rutils::etfenv$VTI
class(pricev); head(pricev)
pricev <- as.data.frame(pricev)
class(pricev); head(pricev)
# Coerce data frame to a data table
data.table::setDT(pricev, keep.rownames=TRUE)
class(pricev); head(pricev)
# Dates are coerced to strings
sapply(pricev, class)
# Coerce xts directly to a data table
dtable <- as.data.table(rutils::etfenv$VTI,
  keep.rownames=TRUE)
class(dtable); head(dtable)
# Dates are not coerced to strings
sapply(dtable, class)
all.equal(pricev, dtable, check.attributes=FALSE)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Package \protect\emph{fst} for High Performance Data Management}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package
      \href{https://cran.r-project.org/web/packages/fst/}{\emph{fst}}
      provides functions for very fast writing and reading of \emph{data frames} from \emph{compressed binary files}.
      \vskip1ex
      The package \emph{fst} writes to \emph{compressed binary files} in the
\texttt{fst} fast-storage format.
      \vskip1ex
      The package \emph{fst} uses the \texttt{LZ4} and \texttt{ZSTD} compression algorithms, and utilizes multithreaded (parallel) processing on multiple CPU cores.
      \vskip1ex
      The package \emph{fst} has extensive documentation:\\
      \hskip1em\url{http://www.fstpackage.org/}
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Install package fst
install.packages("fst")
# Load package fst
library(fst)
# Get documentation for package fst
# Get short description
packageDescription("fst")
# Load help page
help(package="fst")
# List all datasets in "fst"
data(package="fst")
# List all objects in "fst"
ls("package:fst")
# Remove fst from search path
detach("package:fst")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing and Reading Data Using Package \protect\emph{fst}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package
      \href{https://cran.r-project.org/web/packages/fst/}{\emph{fst}}
      allows very fast writing and reading of \emph{data frames} from \emph{compressed binary files} in the \texttt{fst} fast-storage format.
      \vskip1ex
      The function \texttt{fst::write\_fst()} writes to \texttt{.fst} files over \texttt{10} times faster than the function \texttt{write.csv()}, and \texttt{300} times faster than function \texttt{cat()} write to \texttt{.csv} files!
      \vskip1ex
      The function \texttt{fst::fread()} reads from \texttt{.fst} files over \texttt{10} times faster than the function \texttt{read.csv()} from \texttt{.csv} files!
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Read a data frame from CSV file
dirn <- "/Users/jerzy/Develop/lecture_slides/data/"
filen <- file.path(dirn, "weather_delays14.csv")
data.table::setDF(dframe)
class(dframe); dim(dframe)
# Write data frame to .fst file in different ways
fst::write_fst(dframe, path="dframe.fst")
write.csv(dframe, file="dframe2.csv")
# microbenchmark speed of fst::write_fst()
library(microbenchmark)
summary(microbenchmark(
  fst=fst::write_fst(dframe, path="dframe.csv"),
  write_csv=write.csv(dframe, file="dframe2.csv"),
  cat=cat(unlist(dframe), file="dframe3.csv"),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
# fst::read_fst() reads the same data as read.csv()
all.equal(read.csv(filen),
          fst::read_fst("dframe.fst"))
# fst::read_fst() is 10 times faster than read.csv()
summary(microbenchmark(
  fst=fst::read_fst("dframe.fst"),
  read_csv=read.csv(filen),
  times=10))[, c(1, 4, 5)]  # end microbenchmark summary
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Random Access to Large Data Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package
      \href{https://cran.r-project.org/web/packages/fst/}{\emph{fst}}
      allows \emph{random access} to very large \emph{data frames} stored in compressed data files in the \texttt{.fst} format.
      \vskip1ex
      Data frames can be accessed \emph{randomly} by loading only the selected rows and columns into memory, without fully loading the whole data frame.
      \vskip1ex
      function \texttt{fst::fst()} reads an \texttt{.fst} file and returns an \emph{fst\_table} reference object (pointer) to the data, without loading the whole data into memory.
      \vskip1ex
      The \emph{fst\_table} reference provides access to the data similar to a regular \emph{data frame}, but it requires only a small amount of memory because the data isn't loaded into memory.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Coerce TAQ xts to a data frame
library(HighFreq)
taq <- HighFreq::SPY_TAQ
taq <- as.data.frame(taq)
class(taq)
# Coerce data frame to a data table
data.table::setDT(taq, keep.rownames=TRUE)
class(taq); head(taq)
# Get memory size of data table
format(object.size(taq), units="MB")
# Save data table to .fst file
fst::write_fst(taq, path="/Users/jerzy/Develop/data/taq.fst")
# Create reference to .fst file similar to a data frame
refst <- fst::fst("/Users/jerzy/Develop/data/taq.fst")
class(refst)
# Memory size of reference to .fst is very small
format(object.size(refst), units="MB")
# Get sizes of all objects in workspace
sort(sapply(mget(ls()), object.size))
# Reference to .fst can be treated similar to a data table
dim(taq); dim(refst)
fst:::print.fst_table(refst)
# Subset reference to .fst just like a data table
refst[1e4:(1e4+5), ]
      @
  \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Downloading Time Series Data}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{ts} Time Series Using \protect\emph{tseries}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \texttt{get.hist.quote()} can download daily historical data in the \emph{ts} format using the argument \texttt{"retclass="ts"}.
      \vskip1ex
      The default market data provider is \emph{Yahoo} (\texttt{provider="yahoo"}), but \emph{Yahoo} has stopped providing free market data. 
      \vskip1ex
      \texttt{get.hist.quote()} returns a \emph{ts} object with a \texttt{frequency=1}, implying a \emph{"day"} time unit, instead of a \emph{"year"} time unit suitable for \emph{year-fraction} dates.
      \vskip1ex
      The \emph{ts} contains \texttt{NA} values for weekends and holidays.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
load(file="/Users/jerzy/Develop/lecture_slides/data/zoo_data.RData")
library(tseries)  # Load package tseries
# Download MSFT data in ts format
pricev <- suppressWarnings(
  get.hist.quote(
    instrument="MSFT",
    start=Sys.Date()-3*365,
    end=Sys.Date(),
    retclass="ts",
    quote=c("Open","High","Low","Close",
            "AdjClose","Volume"),
    origin="1970-01-01")
)  # end suppressWarnings
# Calculate price adjustment vector
ratio <- as.numeric(pricev[, "AdjClose"]/pricev[, "Close"])
# Adjust OHLC prices
pricadj <- pricev
pricadj[, c("Open","High","Low","Close")] <-
  ratio*pricev[, c("Open","High","Low","Close")]
# Inspect the data
tsp(pricadj)  # frequency=1
head(time(pricadj))
head(pricadj)
tail(pricadj)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{zoo} Time Series Using \protect\emph{tseries}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{get.hist.quote()} downloads historical data from online sources.
      \vskip1ex
      The \texttt{"provider"} argument determines the \emph{online source}, and its default value is \texttt{c("yahoo", "oanda")}.
      \vskip1ex
      The \texttt{"retclass"} argument determines the \emph{return class}, and its default value is \texttt{c("zoo", "its", "ts")}.
      \vskip1ex
      The \texttt{"quote"} argument determines the data fields, and its default value is \texttt{c("Open", "High", "Low", "Close")}.
      \vskip1ex
      The \texttt{"AdjClose"} data field is for the \emph{Close} price adjusted for stock splits and dividends.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(tseries)  # Load package tseries
# Download MSFT data
pricezoo <- suppressWarnings(
  get.hist.quote(
    instrument="MSFT",
    start=Sys.Date()-3*365,
    end=Sys.Date(),
    quote=c("Open","High","Low","Close",
            "AdjClose","Volume"),
    origin="1970-01-01")
)  # end suppressWarnings
      @
      \vspace{-2em}
      <<echo=(-(1:1)),eval=FALSE>>=
load(file="/Users/jerzy/Develop/lecture_slides/data/zoo_data.RData")
class(pricezoo)
dim(pricezoo)
head(pricezoo, 4)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Adjusting \protect\emph{OHLC} Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Stock prices experience jumps due to stock splits and dividends.
      \vskip1ex
      \emph{Adjusted} stock prices are stock prices that have been adjusted so they don't have jumps.
      \vskip1ex
      \emph{OHLC} data can be adjusted for stock splits and dividends.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:2)),eval=FALSE>>=
library(tseries)  # Load package tseries
load(file="/Users/jerzy/Develop/lecture_slides/data/zoo_data.RData")
# Calculate price adjustment vector
ratio <- as.numeric(pricezoo[, "AdjClose"]/pricezoo[, "Close"])
head(ratio, 5)
tail(ratio, 5)
# Adjust OHLC prices
pricedj <- pricezoo
pricedj[, c("Open","High","Low","Close")] <-
  ratio*pricezoo[, c("Open","High","Low","Close")]
head(pricedj)
tail(pricedj)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Data From \protect\emph{Oanda} Using \protect\emph{tseries}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Oanda} is a foreign exchange broker that also provides free historical currency rates data.
      \vskip1ex
      The function \texttt{get.hist.quote()} downloads historical data from online sources.
      \vskip1ex
      The \texttt{"provider"} argument determines the \emph{online source}, and its default value is \texttt{c("yahoo", "oanda")}.
      \vskip1ex
      The \texttt{"retclass"} argument determines the \emph{return class}, and its default value is \texttt{c("zoo", "its", "ts")}.
      \vskip1ex
      The \texttt{"quote"} argument determines the data fields, and its default value is \texttt{c("Open", "High", "Low", "Close")}.
      \vskip1ex
      The function \texttt{complete.cases()} returns \texttt{TRUE} if a row has no \texttt{NA} values.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(tseries)  # Load package tseries
# Download EUR/USD data
priceur <- suppressWarnings(
  get.hist.quote(
    instrument="EUR/USD",
    provider="oanda",
    start=Sys.Date()-3*365,
    end=Sys.Date(),
    origin="1970-01-01")
)  # end suppressWarnings
# Bind and scrub data
pricecombo <- cbind(priceur, pricezoo[, "AdjClose"])
colnames(pricecombo) <- c("EURUSD", "MSFT")
pricecombo <- pricecombo[complete.cases(pricecombo),]
save(pricezoo, pricedj,
     pricev, pricadj,
     priceur, pricecombo,
     file="/Users/jerzy/Develop/lecture_slides/data/zoo_data.RData")
      @
      \vspace{-2em}
      <<echo=(-(1:1)),eval=FALSE>>=
load(file="/Users/jerzy/Develop/lecture_slides/data/zoo_data.RData")
# Inspect the data
class(priceur)
head(priceur, 4)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Stock Prices Using \protect\emph{tseries}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Data for multiple symbols can be downloaded in an \texttt{lapply()} loop, which calls the function \texttt{tseries::get.hist.quote}.
      \vskip1ex
      If the body of an \texttt{apply()} loop returns a \emph{zoo} or \emph{xts} series, then the loop will produce an error, because \texttt{apply()} attempts to coerce its output into a vector or matrix.
      \vskip1ex
      So \texttt{lapply()} should be used instead of \texttt{apply()}.
      \vskip1ex
      The functional \texttt{lapply()} applies a function to a list of objects and returns a list of objects.
      \vskip1ex
      The list of \emph{zoo} time series can be flattened into a single \emph{zoo} series using functions \texttt{do.call()} and \texttt{cbind()}.
      \vskip1ex
      The function \texttt{do.call()} executes a function call using a function name and a list of arguments.
      \vskip1ex
      The function \texttt{do\_call()} from package \emph{rutils} performs the same operation as \texttt{do.call()}, but using recursion, which is much faster and uses less memory.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(tseries)  # Load package tseries
# Download price and volume data for symbolv into list of zoo objects
pricev <- suppressWarnings(
  lapply(symbolv, # Loop for loading data
         get.hist.quote,
         quote=c("AdjClose", "Volume"),
         start=Sys.Date()-3650,
         end=Sys.Date(),
         origin="1970-01-01")  # end lapply
)  # end suppressWarnings
# Flatten list of zoo objects into a single zoo object
pricev <- rutils::do_call(cbind, pricev)
# Or
# pricev <- do.call(cbind, pricev)
# Assign names in format "symboln.Close", "symboln.Volume"
names(pricev) <- as.numeric(sapply(symbolv,
    paste, c("Close", "Volume"), sep="."))
# Save pricev to a comma-separated CSV file
write.zoo(pricev, file="zooseries.csv", sep=",")
# Save pricev to a binary .RData file
save(pricev, file="pricev.RData")
      @
  \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Market Databases}


%%%%%%%%%%%%%%%
\subsection{The \protect\emph{ETF} Database}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Exchange-traded Funds (\emph{ETFs}) are funds which invest in portfolios of assets, such as stocks, commodities, or bonds.
      \vskip1ex
      \emph{ETFs} are shares in portfolios of assets, and they are traded just like stocks.
      \vskip1ex
      \emph{ETFs} provide investors with convenient, low cost, and liquid instruments to invest in various portfolios of assets.
      \vskip1ex
      The file \texttt{etf\_list.csv} contains a database of exchange-traded funds (\emph{ETFs}) and exchange traded notes (\emph{ETNs}).
      \vskip1ex
      We will select a portfolio of \emph{ETFs} for illustrating various investment strategies.
    \column{0.5\textwidth}
    \vspace{-1em}
      <<echo=TRUE,eval=TRUE,size="tiny">>=
# Select ETF symbols for asset allocation
symbolv <- c("SPY", "VTI", "QQQ", "VEU", "EEM", "XLY", "XLP", 
"XLE", "XLF", "XLV", "XLI", "XLB", "XLK", "XLU", "VYM", "IVW", 
"IWB", "IWD", "IWF", "IEF", "TLT", "VNQ", "DBC", "GLD", "USO", 
"VXX", "SVXY", "MTUM", "IVE", "VLUE", "QUAL", "VTV", "USMV", "AIEQ")
# Read etf database into data frame
etflist <- read.csv(file="/Users/jerzy/Develop/lecture_slides/data/etf_list.csv")
rownames(etflist) <- etflist$Symbol
# Select from etflist only those ETF's in symbolv
etflist <- etflist[symbolv, ]
# Shorten names
etfnames <- sapply(etflist$Name, function(name) {
  namesplit <- strsplit(name, split=" ")[[1]]
  namesplit <- namesplit[c(-1, -NROW(namesplit))]
  name_match <- match("Select", namesplit)
  if (!is.na(name_match))
    namesplit <- namesplit[-name_match]
  paste(namesplit, collapse=" ")
})  # end sapply
etflist$Name <- etfnames
etflist["IEF", "Name"] <- "10 year Treasury Bond Fund"
etflist["TLT", "Name"] <- "20 plus year Treasury Bond Fund"
etflist["XLY", "Name"] <- "Consumer Discr. Sector Fund"
etflist["EEM", "Name"] <- "Emerging Market Stock Fund"
etflist["MTUM", "Name"] <- "Momentum Factor Fund"
etflist["SVXY", "Name"] <- "Short VIX Futures"
etflist["VXX", "Name"] <- "Long VIX Futures"
etflist["DBC", "Name"] <- "Commodity Futures Fund"
etflist["USO", "Name"] <- "WTI Oil Futures Fund"
etflist["GLD", "Name"] <- "Physical Gold Fund"
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{\protect\emph{ETF} Database for Investment Strategies}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The database contains \emph{ETFs} representing different \emph{industry sectors} and \emph{investment styles}.
      \vskip1ex
      The \emph{ETFs} with names \emph{X*} represent industry \emph{sector funds} (energy, financial, etc.)
      \vskip1ex
      The \emph{ETFs} with names \emph{I*} represent \emph{style funds} (value, growth, size).
      \vskip1ex
      \emph{IWB} is the Russell 1000 small-cap fund.
      \vskip1ex
      The 
      \href{https://www.ssga.com/us/en/intermediary/etfs/funds/spdr-sp-500-etf-trust-spy}{\emph{SPY ETF}}
      owns the \emph{S\&P500} index constituents.
      \emph{SPY} is the biggest, the most liquid, and the oldest ETF. 
      SPY has over \texttt{\$400} billion of shares outstanding, and trades over \texttt{\$20} billion per day, at a bid-ask spread of only one tick (cent=\texttt{\$0.01}, or about \texttt{0.0022\%}).
      \vskip1ex
      The 
      \href{https://www.invesco.com/qqq-etf/en/home.html}{\emph{QQQ ETF}}
      owns the \emph{Nasdaq-100} index constituents.
      \vskip1ex
      \emph{MTUM} is an \emph{ETF} which owns a stock portfolio representing the \emph{momentum factor}.
      \vskip1ex
      \emph{DBC} is an \emph{ETF} providing the total return on a portfolio of commodity futures.
    \column{0.5\textwidth}
    \vspace{-1em}
      <<echo=FALSE,eval=TRUE,results='asis'>>=
print(xtable::xtable(etflist), comment=FALSE, size="tiny", include.rownames=FALSE)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Exchange Traded Notes (\protect\emph{ETNs})}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{ETNs} are similar to \emph{ETFs}, with the difference that \emph{ETFs} are shares in a fund which owns the underlying assets, while \emph{ETNs} are notes from issuers which promise payouts according to a formula tied to the underlying asset.
      \vskip1ex
      \emph{ETFs} are similar to mutual funds, while \emph{ETNs} are similar to corporate bonds.
      \vskip1ex
      \emph{ETNs} are technically unsecured corporate debt, but instead of fixed coupons, they promise to provide returns on a market index or futures contract.
      \vskip1ex
      The \emph{ETN} issuer promises the payout and is responsible for tracking the index.
      \vskip1ex
      The \emph{ETN} investor has counterparty credit risk to the \emph{ETN} issuer.
    \column{0.5\textwidth}
      \emph{VXX} is an \emph{ETN} providing the total return of \emph{long VIX} futures contracts (specifically the \emph{S\&P} VIX Short-Term Futures Index).
      \vskip1ex
      \emph{VXX} is \emph{bearish} because it's \emph{long} VIX futures, and the VIX \emph{rises} when stock prices \emph{drop}.
      \vskip1ex
      \emph{SVXY} is an \emph{ETF} providing the total return of \emph{short VIX} futures contracts.
      \vskip1ex
      \emph{SVXY} is \emph{bullish} because it's \emph{short} VIX futures, and the VIX \emph{drops} when stock prices \emph{rise}.
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading ETF Prices Using Package \protect\emph{quantmod}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{getSymbols()} downloads time series data into the specified \emph{environment}.
      \vskip1ex
      \texttt{getSymbols()} downloads the daily \emph{OHLC} prices and trading volume (Open, High, Low, Close, Adjusted, Volume).
      \vskip1ex
      \texttt{getSymbols()} creates objects in the specified \emph{environment} from the input strings (names), and assigns the data to those objects, without returning them as a function value, as a \emph{side effect}.
      \vskip1ex
      If the argument \texttt{"auto.assign"} is set to \texttt{FALSE}, then \texttt{getSymbols()} returns the data, instead of assigning it silently.
      \vskip1ex
      \emph{Yahoo} data quality deteriorated significantly in \texttt{2017}, and \emph{Google} data quality is also poor, leaving \emph{Tiingo} and \emph{Alpha Vantage} as the only major providers of free daily \emph{OHLC} stock prices.
      \vskip1ex
      But \href{https://www.quandl.com/}{Quandl} doesn't provide free \emph{ETF} prices, leaving \emph{Alpha Vantage} as the best provider of free daily \emph{ETF} prices.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Select ETF symbols for asset allocation
symbolv <- c("SPY", "VTI", "QQQ", "VEU", "EEM", "XLY", "XLP", 
"XLE", "XLF", "XLV", "XLI", "XLB", "XLK", "XLU", "VYM", "IVW", 
"IWB", "IWD", "IWF", "IEF", "TLT", "VNQ", "DBC", "GLD", "USO", 
"VXX", "SVXY", "MTUM", "IVE", "VLUE", "QUAL", "VTV", "USMV", "AIEQ")
library(rutils)  # Load package rutils
etfenv <- new.env()  # New environment for data
# Boolean vector of symbols already downloaded
isdown <- symbolv %in% ls(etfenv)
# Download data for symbolv using single command - creates pacing error
getSymbols.av(symbolv, adjust=TRUE, env=etfenv,
  output.size="full", api.key="T7JPW54ES8G75310")
# Download data from Alpha Vantage using while loop
nattempts <- 0  # number of download attempts
while ((sum(!isdown) > 0) & (nattempts < 10)) {
  # Download data and copy it into environment
  nattempts <- nattempts + 1
  cat("Download attempt = ", nattempts, "\n")
  for (symboln in na.omit(symbolv[!isdown][1:5])) {
    cat("Processing: ", symboln, "\n")
    tryCatch(  # With error handler
      quantmod::getSymbols.av(symboln, adjust=TRUE, env=etfenv, auto.assign=TRUE, output.size="full", api.key="T7JPW54ES8G75310"),
      # Error handler captures error condition
      error=function(msg) {
        print(paste0("Error handler: ", msg))
      },  # end error handler
      finally=print(paste0("Symbol = ", symboln))
    )  # end tryCatch
  }  # end for
  # Update vector of symbols already downloaded
  isdown <- symbolv %in% ls(etfenv)
  cat("Pausing 1 minute to avoid pacing...\n")
  Sys.sleep(65)
}  # end while
# Download all symbolv using single command - creates pacing error
# quantmod::getSymbols.av(symbolv, env=etfenv, adjust=TRUE, from="2005-01-03", output.size="full", api.key="T7NHW54ES8GG501C")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Inspecting ETF Prices in an Environment}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{get()} retrieves objects that are referenced using character strings, instead of their names.
      \vskip1ex
      The function \texttt{eapply()} is similar to \texttt{lapply()}, and applies a function to objects in an \emph{environment}, and returns a list.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
ls(etfenv)  # List files in etfenv
# Get class of object in etfenv
class(get(x=symbolv[1], envir=etfenv))
# Another way
class(etfenv$VTI)
colnames(etfenv$VTI)
# Get first 3 rows of data
head(etfenv$VTI, 3)
# Get last 11 rows of data
tail(etfenv$VTI, 11)
# Get class of all objects in etfenv
eapply(etfenv, class)
# Get class of all objects in R workspace
lapply(ls(), function(namev) class(get(namev)))
# Get end dates of all objects in etfenv
as.Date(sapply(etfenv, end))
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Adjusting Stock Prices Using Package \protect\emph{quantmod}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Traded stock and bond prices experience jumps after splits and dividends, and must be adjusted to account for them.
      \vskip1ex
      The function \texttt{adjustOHLC()} adjusts \emph{OHLC} prices.
      \vskip1ex
      The function \texttt{get()} retrieves objects that are referenced using character strings, instead of their names.
      \vskip1ex
      The function \texttt{assign()} assigns a value to an object in a specified \emph{environment}, by referencing it using a character string (name).
      \vskip1ex
      The functions \texttt{get()} and \texttt{assign()} allow retrieving and assigning values to objects that are referenced using character strings.
      \vskip1ex
      The function \texttt{mget()} accepts a vector of strings and returns a list of the corresponding objects extracted from an \emph{environment}.
      \vskip1ex
      If the argument \texttt{"adjust"} in function \texttt{getSymbols()} is set to \texttt{TRUE}, then \texttt{getSymbols()} returns adjusted data.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
# Check of object is an OHLC time series
is.OHLC(etfenv$VTI)
# Adjust single OHLC object using its name
etfenv$VTI <- adjustOHLC(etfenv$VTI, use.Adjusted=TRUE)

# Adjust OHLC object using string as name
assign(symbolv[1], adjustOHLC(
    get(x=symbolv[1], envir=etfenv), use.Adjusted=TRUE),
  envir=etfenv)

# Adjust objects in environment using vector of strings
for (symboln in ls(etfenv)) {
  assign(symboln,
         adjustOHLC(get(symboln, envir=etfenv), use.Adjusted=TRUE),
         envir=etfenv)
}  # end for
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Extracting Time Series from Environments}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{mget()} accepts a vector of strings and returns a list of the corresponding objects extracted from an \emph{environment}.
      \vskip1ex
      The extractor (accessor) functions from package \emph{quantmod}: \texttt{Cl()}, \texttt{Vo()}, etc., extract columns from \emph{OHLC} data.
      \vskip1ex
      A list of \emph{xts} series can be flattened into a single \emph{xts} series using the function \texttt{do.call()}.
      \vskip1ex
      The function \texttt{do.call()} executes a function call using a function name and a list of arguments.
      \vskip1ex
      \texttt{do.call()} passes the list elements individually, instead of passing the whole list as one argument.
      \vskip1ex
      The function \texttt{eapply()} is similar to \texttt{lapply()}, and applies a function to objects in an \emph{environment}, and returns a list.
      \vskip1ex
      Time series can also be extracted from an \emph{environment} by coercing it into a \texttt{list}, and then subsetting and merging it into an \emph{xts} series using the function \texttt{do.call()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
library(rutils)  # Load package rutils
# Define ETF symbols
symbolv <- c("VTI", "VEU", "IEF", "VNQ")
# Extract symbolv from rutils::etfenv
pricev <- mget(symbolv, envir=rutils::etfenv)
# pricev is a list of xts series
class(pricev)
class(pricev[[1]])
tail(pricev[[1]])
# Extract close prices
pricev <- lapply(pricev, quantmod::Cl)
# Collapse list into time series the hard way
prices2 <- cbind(pricev[[1]], pricev[[2]], pricev[[3]], pricev[[4]])
class(price2)
dim(price2)
# Collapse list into time series using do.call()
pricev <- do.call(cbind, pricev)
all.equal(price2, pricev)
class(pricev)
dim(pricev)
# Or extract and cbind in single step
pricev <- do.call(cbind, lapply(
  mget(symbolv, envir=rutils::etfenv), quantmod::Cl))
# Or extract and bind all data, subset by symbolv
pricev <- lapply(symbolv, function(symboln) {
    quantmod::Cl(get(symboln, envir=rutils::etfenv))
})  # end lapply
# Or loop over etfenv without anonymous function
pricev <- do.call(cbind,
  lapply(as.list(rutils::etfenv)[symbolv], quantmod::Cl))
# Same, but works only for OHLC series - produces error
pricev <- do.call(cbind,
  eapply(rutils::etfenv, quantmod::Cl)[symbolv])
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Managing Time Series}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Time series columns can be renamed, and then saved into \texttt{.csv} files.
      \vskip1ex
      The function \texttt{strsplit()} splits the elements of a character vector.
      \vskip1ex
      The package \emph{zoo} contains functions \texttt{write.zoo()} and \texttt{read.zoo()} for writing and reading \emph{zoo} time series from \texttt{.txt} and \texttt{.csv} files.
      \vskip1ex
      The function \texttt{eapply()} is similar to \texttt{lapply()}, and applies a function to objects in an \emph{environment}, and returns a list.
      \vskip1ex
      The function \texttt{assign()} assigns a value to an object in a specified \emph{environment}, by referencing it using a character string (name).
      \vskip1ex
      The function \texttt{save()} writes objects to compressed binary \texttt{.RData} files.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Column names end with ".Close"
colnames(pricev)
strsplit(colnames(pricev), split="[.]")
do.call(rbind, strsplit(colnames(pricev), split="[.]"))
do.call(rbind, strsplit(colnames(pricev), split="[.]"))[, 1]
# Drop ".Close" from colnames
colnames(pricev) <- rutils::get_name(colnames(pricev))
# Or
# colnames(pricev) <- do.call(rbind,
#   strsplit(colnames(pricev), split="[.]"))[, 1]
tail(pricev, 3)
# Which objects in global environment are class xts?
unlist(eapply(globalenv(), is.xts))
# Save xts to csv file
write.zoo(pricev,
  file="/Users/jerzy/Develop/lecture_slides/data/etf_series.csv", sep=",")
# Copy prices into etfenv
etfenv$prices <- pricev
# Or
assign("pricev", pricev, envir=etfenv)
# Save to .RData file
save(etfenv, file="etf_data.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Calculating Percentage Returns from Close Prices}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{quantmod::dailyReturn()} calculates the percentage daily returns from the \emph{Close} prices.
      \vskip1ex
      The \texttt{lapply()} and \texttt{sapply()} functionals perform a loop over the columns of \emph{zoo} and \emph{xts} series.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Extract VTI prices
pricev <- etfenv$prices[ ,"VTI"]
pricev <- na.omit(pricev)
# Calculate percentage returns "by hand"
pricel <- as.numeric(pricev)
pricel <- c(pricel[1], pricel[-NROW(pricel)])
pricel <- xts(pricel, zoo::index(pricev))
retp <- (pricev-pricel)/pricel
# Calculate percentage returns using dailyReturn()
retd <- quantmod::dailyReturn(pricev)
head(cbind(retd, retp))
all.equal(retd, retp, check.attributes=FALSE)
# Calculate returns for all prices in etfenv$prices
retp <- lapply(etfenv$prices, function(xtsv) {
  retd <- quantmod::dailyReturn(na.omit(xtsv))
  colnames(retd) <- names(xtsv)
  retd
})  # end lapply
# "retp" is a list of xts
class(retp)
class(retp[[1]])
# Flatten list of xts into a single xts
retp <- do.call(cbind, retp)
class(retp)
dim(retp)
# Copy retp into etfenv and save to .RData file
# assign("retp", retp, envir=etfenv)
etfenv$retp <- retp
save(etfenv, file="/Users/jerzy/Develop/lecture_slides/data/etf_data.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Managing Data Inside Environments}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{as.environment()} coerces objects (listv) into an environment.
      \vskip1ex
      The function \texttt{eapply()} is similar to \texttt{lapply()}, and applies a function to objects in an \emph{environment}, and returns a list.
      \vskip1ex
      The function \texttt{mget()} accepts a vector of strings and returns a list of the corresponding objects extracted from an \emph{environment}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
library(rutils)
startd <- "2012-05-10"; endd <- "2013-11-20"
# Select all objects in environment and return as environment
newenv <- as.environment(eapply(etfenv, "[",
                  paste(startd, endd, sep="/")))
# Select only symbolv in environment and return as environment
newenv <- as.environment(
  lapply(as.list(etfenv)[symbolv], "[",
         paste(startd, endd, sep="/")))
# Extract and cbind Close prices and return to environment
assign("prices", rutils::do_call(cbind,
  lapply(ls(etfenv), function(symboln) {
    xtsv <- quantmod::Cl(get(symboln, etfenv))
    colnames(xtsv) <- symboln
    xtsv
  })), envir=newenv)
# Get sizes of OHLC xts series in etfenv
sapply(mget(symbolv, envir=etfenv), object.size)
# Extract and cbind adjusted prices and return to environment
colname <- function(xtsv)
  strsplit(colnames(xtsv), split="[.]")[[1]][1]
assign("prices", rutils::do_call(cbind,
               lapply(mget(etfenv$symbolv, envir=etfenv),
                      function(xtsv) {
                        xtsv <- Ad(xtsv)
                        colnames(xtsv) <- colname(xtsv)
                        xtsv
               })), envir=newenv)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Stock Databases And Survivorship Bias}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The file \texttt{sp500\_constituents.csv} contains a \emph{data frame} of over \texttt{700} present (and also some past) \emph{S\&P500} index constituents.
      \vskip1ex
      The file \texttt{sp500\_constituents.csv} is updated with stocks recently added to the \emph{S\&P500} index by downloading the
      \href{https://www.ssga.com/us/en/intermediary/etfs/funds/spdr-sp-500-etf-trust-spy}{\emph{SPY ETF Holdings}}.
      \vskip1ex
      But the file \texttt{sp500\_constituents.csv} doesn't include companies that have gone bankrupt.
      For example, it doesn't include Enron, which was in the \emph{S\&P500} index before it went bankrupt in 2001.
      \vskip1ex
      Most databases of stock prices don't include companies that have gone bankrupt or have been liquidated.
      \vskip1ex
      This introduces a \emph{survivorship bias} to the data, which can skew portfolio simulations and strategy backtests.
      \vskip1ex
      Accurate strategy simulations require starting with a portfolio of companies at a "point in time" in the past, and tracking them over time.
      \vskip1ex
      Research databases like the 
      \href{https://wrds-www.wharton.upenn.edu}{\emph{WRDS}} 
      database provide stock prices of companies that are no longer traded.
      \vskip1ex
      The stock tickers are stored in the column \texttt{"Ticker"} of the \texttt{sp500} \emph{data frame}.
      \vskip1ex
      Some tickers (like "BRK.B" and "BF.B") are not valid symbols in \emph{Tiingo}, so they must be renamed.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Load data frame of S&P500 constituents from CSV file
sp500 <- read.csv(file="/Users/jerzy/Develop/lecture_slides/data/sp500_constituents.csv")
# Inspect data frame of S&P500 constituents
dim(sp500)
colnames(sp500)
# Extract tickers from the column Ticker
symbolv <- sp500$Ticker
# Get duplicate tickers
tablev <- table(symbolv)
duplicatv <- tablev[tablev > 1]
duplicatv <- names(duplicatv)
# Get duplicate records (rows) of sp500
sp500[symbolv %in% duplicatv, ]
# Get unique tickers
symbolv <- unique(symbolv)
# Find index of ticker "BRK.B"
which(symbolv=="BRK.B")
# Rename "BRK.B" to "BRK-B" and "BF.B" to "BF-B"
symbolv[which(symbolv=="BRK.B")] <- "BRK-B"
symbolv[which(symbolv=="BF.B")] <- "BF-B"
@
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Stock Time Series From \protect\emph{Tiingo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Yahoo} data quality deteriorated significantly in \texttt{2017}, and \emph{Google} data quality is also poor, leaving \emph{Tiingo}, \emph{Alpha Vantage}, and \href{https://www.quandl.com/}{Quandl} as the only major providers of free daily \emph{OHLC} stock prices.
      \vskip1ex
      But \href{https://www.quandl.com/}{Quandl} doesn't provide free \emph{ETF} prices, while \emph{Tiingo} does.
      \vskip1ex
      The function \texttt{getSymbols()} has a \emph{method} for downloading time series data from \emph{Tiingo}, called \texttt{getSymbols.tiingo()}.
      \vskip1ex
      Users must first obtain a \emph{Tiingo} \emph{API key}, and then pass it in \texttt{getSymbols.tiingo()} calls:\\
      https://www.tiingo.com/
      \vskip1ex
      Note that the data are downloaded as \texttt{xts} time series, with a date-time index of class \texttt{POSIXct} (not \texttt{Date}).
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Load package rutils
library(rutils)
# Create new environment for data
sp500env <- new.env()
# Boolean vector of symbols already downloaded
isdown <- symbolv %in% ls(sp500env)
# Download in while loop from Tiingo and copy into environment
nattempts <- 0  # Number of download attempts
while ((sum(!isdown) > 0) & (nattempts<3)) {
  # Download data and copy it into environment
  nattempts <- nattempts + 1
  cat("Download attempt = ", nattempts, "\n")
  for (symboln in symbolv[!isdown]) {
    cat("processing: ", symboln, "\n")
    tryCatch(  # With error handler
      quantmod::getSymbols(symboln, src="tiingo", adjust=TRUE, auto.assign=TRUE,
                 from="1990-01-01", env=sp500env, api.key="j84ac2b9c5bde2d68e33034f65d838092c6c9f10"),
      # Error handler captures error condition
      error=function(msg) {
        print(paste0("Error handler: ", msg))
      },  # end error handler
      finally=print(paste0("Symbol = ", symboln))
    )  # end tryCatch
  }  # end for
  # Update vector of symbols already downloaded
  isdown <- symbolv %in% ls(sp500env)
  Sys.sleep(2)  # Wait 2 seconds until next attempt
}  # end while
class(sp500env$AAPL)
class(zoo::index(sp500env$AAPL))
tail(sp500env$AAPL)
symbolv[!isdown]
@
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Coercing Date-time Indices}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The date-time indices of the \emph{OHLC} stock prices are in the \texttt{POSIXct} format suitable for intraday prices, not daily prices.
      \vskip1ex
      The function \texttt{as.Date()} coerces \texttt{POSIXct} objects into \texttt{Date} objects.
      \vskip1ex
      The function \texttt{get()} retrieves objects that are referenced using character strings, instead of their names.
      \vskip1ex
      The function \texttt{assign()} assigns a value to an object in a specified \emph{environment}, by referencing it using a character string (name).
      \vskip1ex
      The functions \texttt{get()} and \texttt{assign()} allow retrieving and assigning values to objects that are referenced using character strings.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# The date-time index of AAPL is POSIXct
class(zoo::index(sp500env$AAPL))
# Coerce the date-time index of AAPL to Date
zoo::index(sp500env$AAPL) <- as.Date(zoo::index(sp500env$AAPL))
# Coerce all the date-time indices to Date
for (symboln in ls(sp500env)) {
  ohlc <- get(symboln, envir=sp500env)
  zoo::index(ohlc) <- as.Date(zoo::index(ohlc))
  assign(symboln, ohlc, envir=sp500env)
}  # end for
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Managing Exceptions in Stock Symbols}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The column names for symbol \texttt{"LOW"} (Lowe's company) must be renamed for the extractor function \texttt{quantmod::Lo()} to work properly.
      \vskip1ex
      Tickers which contain a dot in their name (like "BRK.B") are not valid symbols in \texttt{R}, so they must be downloaded separately and renamed.
      <<echo=TRUE,eval=FALSE>>=
# "LOW.Low" is a bad column name
colnames(sp500env$LOW)
strsplit(colnames(sp500env$LOW), split="[.]")
do.call(cbind, strsplit(colnames(sp500env$LOW), split="[.]"))
do.call(cbind, strsplit(colnames(sp500env$LOW), split="[.]"))[2, ]
# Extract proper names from column names
namev <- rutils::get_name(colnames(sp500env$LOW), field=2)
# Or
# namev <- do.call(rbind, strsplit(colnames(sp500env$LOW),
#                                   split="[.]"))[, 2]
# Rename "LOW" colnames to "LOWES"
colnames(sp500env$LOW) <- paste("LOVES", namev, sep=".")
sp500env$LOWES <- sp500env$LOW
rm(LOW, envir=sp500env)
# Rename BF-B colnames to "BFB"
colnames(sp500env$"BF-B") <- paste("BFB", namev, sep=".")
sp500env$BFB <- sp500env$"BF-B"
rm("BF-B", envir=sp500env)
# Rename BRK-B colnames
sp500env$BRKB <- sp500env$`BRK-B`
rm(`BRK-B`, envir=sp500env)
colnames(sp500env$BRKB) <- gsub("BRK-B", "BRKB", colnames(sp500env$BRKB))
# Save OHLC prices to .RData file
save(sp500env, file="/Users/jerzy/Develop/lecture_slides/data/sp500.RData")
# Download "BRK.B" separately with auto.assign=FALSE
# BRKB <- quantmod::getSymbols("BRK-B", auto.assign=FALSE, src="tiingo", adjust=TRUE, from="1990-01-01", api.key="j84ac2b9c5bde2d68e33034f65d838092c6c9f10")
# colnames(BRKB) <- paste("BRKB", namev, sep=".")
# sp500env$BRKB <- BRKB
      @
    \column{0.5\textwidth}
      \includegraphics[width=0.5\paperwidth]{figure/lowes_stock.png}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Plot OHLC candlestick chart for LOWES
chart_Series(x=sp500env$LOWES["2019-12/"],
  TA="add_Vo()", name="LOWES OHLC Stock Prices")
# Plot dygraph
dygraphs::dygraph(sp500env$LOWES["2019-12/", -5], main="LOWES OHLC Stock Prices") %>%
  dyCandlestick()
      @

  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{\protect\emph{S\&P500} Stock Index Constituent Prices}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The file \texttt{sp500.RData} contains the \emph{environment} \texttt{sp500\_env} with \emph{OHLC} prices and trading volumes of \emph{S\&P500} stock index constituents.
      \vskip1ex
      The \emph{S\&P500} stock index constituent data is of poor quality before \texttt{2000}, so we'll mostly use the data after the year \texttt{2000}.
      <<echo=TRUE,eval=FALSE>>=
# Load S&P500 constituent stock prices
load("/Users/jerzy/Develop/lecture_slides/data/sp500.RData")
pricev <- eapply(sp500env, quantmod::Cl)
pricev <- rutils::do_call(cbind, pricev)
# Carry forward non-NA prices
pricev <- zoo::na.locf(pricev, na.rm=FALSE)
# Drop ".Close" from column names
colnames(pricev)
colnames(pricev) <- rutils::get_name(colnames(pricev))
# Or
# colnames(pricev) <- do.call(rbind,
#   strsplit(colnames(pricev), split="[.]"))[, 1]
# Calculate percentage returns of the S&P500 constituent stocks
# retp <- xts::diff.xts(log(pricev))
retp <- xts::diff.xts(pricev)/
  rutils::lagit(pricev, pad_zeros=FALSE)
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
samplev <- sample(NCOL(retp), s=100, replace=FALSE)
prices100 <- pricev[, samplev]
returns100 <- retp[, samplev]
save(pricev, prices100,
  file="/Users/jerzy/Develop/lecture_slides/data/sp500_prices.RData")
save(retp, returns100,
  file="/Users/jerzy/Develop/lecture_slides/data/sp500_returns.RData")
      @
    \column{0.5\textwidth}
      \includegraphics[width=0.5\paperwidth]{figure/sp500_without_prices.png}
      <<echo=TRUE,eval=FALSE>>=
# Calculate number of constituents without prices
datav <- rowSums(is.na(pricev))
datav <- xts::xts(datav, order.by=zoo::index(pricev))
dygraphs::dygraph(datav, main="Number of S&P500 Constituents Without Prices") %>%
  dyOptions(colors="blue", strokeWidth=2)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{\protect\emph{S\&P500} Stock Portfolio Index}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The price-weighted index of \emph{S\&P500} constituents closely follows the VTI \emph{ETF}.
      <<echo=TRUE,eval=FALSE>>=
# Calculate price weighted index of constituent
ncols <- NCOL(pricev)
pricev <- zoo::na.locf(pricev, fromLast=TRUE)
indeks <- xts(rowSums(pricev)/ncols, zoo::index(pricev))
colnames(indeks) <- "index"
# Combine index with VTI
datav <- cbind(indeks[zoo::index(etfenv$VTI)], etfenv$VTI[, 4])
colv <- c("index", "VTI")
colnames(datav) <- colv
# Plot index with VTI
endd <- rutils::calc_endpoints(datav, interval="weeks")
dygraphs::dygraph(log(datav)[endd],
  main="S&P 500 Price-weighted Index and VTI") %>%
  dyAxis("y", label=colv[1], independentTicks=TRUE) %>%
  dyAxis("y2", label=colv[2], independentTicks=TRUE) %>%
  dySeries(name=colv[1], axis="y", col="red") %>%
  dySeries(name=colv[2], axis="y2", col="blue")
      @
    \column{0.5\textwidth}
      \includegraphics[width=0.5\paperwidth]{figure/sp500_portfolio_index.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Writing Time Series To Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The data from \emph{Tiingo} is downloaded as \texttt{xts} time series, with a date-time index of class \texttt{POSIXct} (not \texttt{Date}).
      \vskip1ex
      The function \texttt{save()} writes objects to compressed binary \texttt{.RData} files.
      \vskip1ex
      The easiest way to share data between \texttt{R} and \texttt{Excel} is through \texttt{.csv} files.
      \vskip1ex
      The package \emph{zoo} contains functions \texttt{write.zoo()} and \texttt{read.zoo()} for writing and reading \emph{zoo} time series from \texttt{.txt} and \texttt{.csv} files.
      \vskip1ex
      The function \texttt{data.table::fread()} reads from \texttt{.csv} files over \texttt{6} times faster than the function \texttt{read.csv()}!
      \vskip1ex
      The function \texttt{data.table::fwrite()} writes to \texttt{.csv} files over \texttt{12} times faster than the function \texttt{write.csv()}, and \texttt{278} times faster than function \texttt{cat()}!
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Save the environment to compressed .RData file
dirn <- "/Users/jerzy/Develop/lecture_slides/data/"
save(sp500env, file=paste0(dirn, "sp500.RData"))
# Save the ETF prices into CSV files
dirn <- "/Users/jerzy/Develop/lecture_slides/data/SP500/"
for (symboln in ls(sp500env)) {
  zoo::write.zoo(sp500env$symbol, file=paste0(dirn, symboln, ".csv"))
}  # end for
# Or using lapply()
filens <- lapply(ls(sp500env), function(symboln) {
  xtsv <- get(symboln, envir=sp500env)
  zoo::write.zoo(xtsv, file=paste0(dirn, symboln, ".csv"))
  symboln
})  # end lapply
unlist(filens)
# Or using eapply() and data.table::fwrite()
filens <- eapply(sp500env , function(xtsv) {
  filen <- rutils::get_name(colnames(xtsv)[1])
  data.table::fwrite(data.table::as.data.table(xtsv), file=paste0(dirn, filen, ".csv"))
  filen
})  # end eapply
unlist(filens)
@
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Reading Time Series from Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{load()} reads data from \texttt{.RData} files, and \emph{invisibly} returns a vector of names of objects created in the workspace.
      \vskip1ex
      The function \texttt{Sys.glob()} listv files matching names obtained from wildcard expansion.
      \vskip1ex
      The easiest way to share data between \texttt{R} and \texttt{Excel} is through \texttt{.csv} files.
      \vskip1ex
      The function \texttt{as.Date()} parses \texttt{character} strings, and coerces \texttt{numeric} and \texttt{POSIXct} objects into \texttt{Date} objects.
      \vskip1ex
      The function \texttt{data.table::setDF()} coerces a \emph{data table} object into a \emph{data frame} using a \emph{side effect}, without making copies of data.
      \vskip1ex
      The function \texttt{data.table::fread()} reads from \texttt{.csv} files over \texttt{6} times faster than the function \texttt{read.csv()}!
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Load the environment from compressed .RData file
dirn <- "/Users/jerzy/Develop/lecture_slides/data/"
load(file=paste0(dirn, "sp500.RData"))
# Get all the .csv file names in the directory
dirn <- "/Users/jerzy/Develop/lecture_slides/data/SP500/"
filens <- Sys.glob(paste0(dirn, "*.csv"))
# Create new environment for data
sp500env <- new.env()
for (filen in filens) {
  xtsv <- xts::as.xts(zoo::read.csv.zoo(filen))
  symboln <- rutils::get_name(colnames(xtsv)[1])
  # symboln <- strsplit(colnames(xtsv), split="[.]")[[1]][1]
  assign(symboln, xtsv, envir=sp500env)
}  # end for
# Or using fread()
for (filen in filens) {
  xtsv <- data.table::fread(filen)
  data.table::setDF(xtsv)
  xtsv <- xts::xts(xtsv[, -1], as.Date(xtsv[, 1]))
  symboln <- rutils::get_name(colnames(xtsv)[1])
  assign(symboln, xtsv, envir=sp500env)
}  # end for
@
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Stock Time Series From \protect\emph{Alpha Vantage}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Yahoo} data quality deteriorated significantly in \texttt{2017}, and \emph{Google} data quality is also poor, leaving \emph{Tiingo}, \emph{Alpha Vantage}, and \href{https://www.quandl.com/}{Quandl} as the only major providers of free daily \emph{OHLC} stock prices.
      \vskip1ex
      But \href{https://www.quandl.com/}{Quandl} doesn't provide free \emph{ETF} prices, while \emph{Alpha Vantage} does.
      \vskip1ex
      The function \texttt{getSymbols()} has a \emph{method} for downloading time series data from \emph{Alpha Vantage}, called \texttt{getSymbols.av()}.
      \vskip1ex
      Users must first obtain an \emph{Alpha Vantage} \emph{API key}, and then pass it in \texttt{getSymbols.av()} calls:\\
      https://www.alphavantage.co/
      \vskip1ex
      The function \texttt{adjustOHLC()} with argument \texttt{use.Adjusted=TRUE}, adjusts all the \emph{OHLC} price columns, using the \emph{Adjusted} price column.
    \column{0.5\textwidth}
      \vspace{-2em}
      <<echo=TRUE,eval=FALSE>>=
# Remove all files from environment(if necessary)
rm(list=ls(sp500env), envir=sp500env)
# Download in while loop from Alpha Vantage and copy into environment
isdown <- symbolv %in% ls(sp500env)
nattempts <- 0
while ((sum(!isdown) > 0) & (nattempts < 10)) {
  # Download data and copy it into environment
  nattempts <- nattempts + 1
  for (symboln in symbolv[!isdown]) {
    cat("processing: ", symboln, "\n")
    tryCatch(  # With error handler
      quantmod::getSymbols(symboln, src="av", adjust=TRUE, auto.assign=TRUE, env=sp500env,
                 output.size="full", api.key="T7JPW54ES8G75310"),
      # error handler captures error condition
      error=function(msg) {
        print(paste0("Error handler: ", msg))
      },  # end error handler
      finally=print(paste0("Symbol = ", symboln))
    )  # end tryCatch
  }  # end for
  # Update vector of symbols already downloaded
  isdown <- symbolv %in% ls(sp500env)
  Sys.sleep(2)  # Wait 2 seconds until next attempt
}  # end while
# Adjust all OHLC prices in environment
for (symboln in ls(sp500env)) {
  assign(symboln,
    adjustOHLC(get(x=symboln, envir=sp500env), use.Adjusted=TRUE),
    envir=sp500env)
}  # end for
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading The \protect\emph{S\&P500} Index Time Series From \protect\emph{Yahoo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{S\&P500} stock market index is a capitalization-weighted average of the 500 largest U.S. companies, and covers about 80\% of the U.S. stock market capitalization.
      \vskip1ex
      Notice: 
      \href{https://algotrading101.com/learn/yahoo-finance-api-guide/}{\emph{Yahoo} no longer provides a public API for data}.
      \vskip1ex
      There are workarounds but they're tedious.
      \vskip1ex
      \emph{Yahoo} provides daily \emph{OHLC} prices for the \emph{S\&P500} index (symbol \emph{\textasciicircum{}GSPC}), and for the \emph{S\&P500} total return index (symbol \emph{\textasciicircum{}SP500TR}).
      \vskip1ex
      But special characters in some stock symbols, like \texttt{"-"} or \texttt{"\textasciicircum{}"} are not allowed in \texttt{R} names.
      \vskip1ex
      For example, the symbol \emph{\textasciicircum{}GSPC} for the \emph{S\&P500} stock market index isn't a valid name in \texttt{R}.
      \vskip1ex
      The function \texttt{setSymbolLookup()} creates valid names corresponding to stock symbols, which are then used by the function \texttt{getSymbols()} to create objects with the valid names.
      \vskip1ex
      \emph{Yahoo} data quality deteriorated significantly in \texttt{2017}, and \emph{Google} data quality is also poor, leaving \emph{Alpha Vantage} and \href{https://www.quandl.com/}{Quandl} as the only major providers of free daily \emph{OHLC} stock prices.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
# Assign name SP500 to ^GSPC symbol
quantmod::setSymbolLookup(SP500=list(name="^GSPC", src="yahoo"))
quantmod::getSymbolLookup()
# View and clear options
options("getSymbols.sources")
options(getSymbols.sources=NULL)
# Download S&P500 prices into etfenv
quantmod::getSymbols("SP500", env=etfenv,
    adjust=TRUE, auto.assign=TRUE, from="1990-01-01")

chart_Series(x=etfenv$SP500["2016/"],
             TA="add_Vo()", name="S&P500 index")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading The \protect\emph{DJIA} Index Time Series From \protect\emph{Yahoo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The Dow Jones Industrial Average (\emph{DJIA}) stock market index is a price-weighted average of the 30 largest U.S. companies (same number of shares per company).
      \vskip1ex
      \emph{Yahoo} provides daily \emph{OHLC} prices for the \emph{DJIA} index (symbol \emph{\textasciicircum{}DJI}), and for the \emph{DJITR} total return index (symbol \emph{DJITR}).
      \vskip1ex
      But special characters in some stock symbols, like \texttt{"-"} or \texttt{"\textasciicircum{}"} are not allowed in \texttt{R} names.
      \vskip1ex
      For example, the symbol \emph{\textasciicircum{}DJI} for the \emph{DJIA} stock market index isn't a valid name in \texttt{R}.
      \vskip1ex
      The function \texttt{setSymbolLookup()} creates valid names corresponding to stock symbols, which are then used by the function \texttt{getSymbols()} to create objects with the valid names.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
# Assign name DJIA to ^DJI symbol
setSymbolLookup(DJIA=list(name="^DJI", src="yahoo"))
getSymbolLookup()
# view and clear options
options("getSymbols.sources")
options(getSymbols.sources=NULL)
# Download DJIA prices into etfenv
quantmod::getSymbols("DJIA", env=etfenv,
    adjust=TRUE, auto.assign=TRUE, from="1990-01-01")
chart_Series(x=etfenv$DJIA["2016/"],
             TA="add_Vo()", name="DJIA index")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Calculating Prices and Returns From \protect\emph{OHLC} Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{na.locf()} from package \emph{zoo} replaces \texttt{NA} values with the most recent non-\texttt{NA} values prior to it.
      \vskip1ex
      The function \texttt{na.locf()} with argument \texttt{fromLast=TRUE} replaces \texttt{NA} values with non-\texttt{NA} values in reverse order, starting from the end.
      \vskip1ex
      The function \texttt{rutils::get\_name()} extracts symbol names (tickers) from a vector of character strings.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
# Calculate prices from OHLC data of the S&P500 stocks
pricev <- eapply(sp500env, quantmod::Cl)
pricev <- rutils::do_call(cbind, pricev)
# Carry forward non-NA prices
pricev <- zoo::na.locf(pricev, na.rm=FALSE)
# Get first column name
colnames(pricev[, 1])
rutils::get_name(colnames(pricev[, 1]))
# Modify column names
colnames(pricev) <- rutils::get_name(colnames(pricev))
# Or
# colnames(pricev) <- do.call(rbind,
#   strsplit(colnames(pricev), split="[.]"))[, 1]
# Calculate percentage returns
retp <- xts::diff.xts(pricev)/
  rutils::lagit(pricev, pad_zeros=FALSE)
# Select a random sample of 100 prices and returns
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
samplev <- sample(NCOL(retp), s=100, replace=FALSE)
prices100 <- pricev[, samplev]
returns100 <- retp[, samplev]
# Save the data into binary files
save(pricev, prices100,
     file="/Users/jerzy/Develop/lecture_slides/data/sp500_prices.RData")
save(retp, returns100,
     file="/Users/jerzy/Develop/lecture_slides/data/sp500_returns.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Stock Prices From Polygon}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \href{http://polygon.io/}{\emph{Polygon}}
      is a premium provider of live and historical stock price data, both daily and intraday (minutes). 
      \vskip1ex
      \emph{Polygon} provides $2$ years of daily historical stock prices for free. 
      But users must first obtain a Polygon \emph{API key}.
      \vskip1ex
      \emph{Polygon} provides the historical \emph{OHLC} stock prices in \emph{JSON} format.
      \vskip1ex
      \emph{JSON} (JavaScript Object Notation) is a data format consisting of symbol-value pairs.
      \vskip1ex
      The package
      \href{https://cran.r-project.org/web/packages/jsonlite/}{\emph{jsonlite}}
      contains functions for managing data in \emph{JSON} format.
      \vskip1ex
      The functions \texttt{fromJSON()} and \texttt{toJSON()} convert data from \emph{JSON} format to \texttt{R} objects, and vice versa.
      \vskip1ex
      The functions \texttt{read\_json()} and \texttt{write\_json()} read and write \emph{JSON} format data in files.
      \vskip1ex
      The function \texttt{download.file()} downloads data from an internet website \texttt{URL} and writes it to a file.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Setup code
symboln <- "SPY"
startd <- as.Date("1990-01-01")
todayd <- Sys.Date()
tspan <- "day"
# Replace below your own Polygon API key
apikey <- "SEpnsBpiRyONMJdl48r6dOo0_pjmCu5r"
# Create url for download
urll <- paste0("https://api.polygon.io/v2/aggs/ticker/", symboln, "/range/1/", tspan, "/", startd, "/", todayd, "?adjusted=true&sort=asc&limit=50000&apiKey=", apikey)
# Download SPY OHLC prices in JSON format from Polygon
ohlc <- jsonlite::read_json(urll)
class(ohlc)
NROW(ohlc)
names(ohlc)
# Extract list of prices from json object
ohlc <- ohlc$results
# Coerce from list to matrix
ohlc <- lapply(ohlc, unlist)
ohlc <- do.call(rbind, ohlc)
# Coerce time from milliseconds to dates
datev <- ohlc[, "t"]/1e3
datev <- as.POSIXct(datev, origin="1970-01-01")
datev <- as.Date(datev)
tail(datev)
# Coerce from matrix to xts
ohlc <- ohlc[, c("o","h","l","c","v","vw")]
colnames(ohlc) <- c("Open", "High", "Low", "Close", "Volume", "VWAP")
ohlc <- xts::xts(ohlc, order.by=datev)
tail(ohlc)
# Save the xts time series to compressed RData file
save(ohlc, file="/Users/jerzy/Data/spy_daily.RData")
# Candlestick plot of SPY OHLC prices
dygraphs::dygraph(ohlc[, 1:4], main=paste("Candlestick Plot of", symboln, "OHLC prices")) %>% 
  dygraphs::dyCandlestick()
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Multiple Stock Prices From Polygon}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.4\textwidth}
      The stock prices for multiple stocks can be downloaded in a \texttt{while()} loop.
      <<echo=TRUE,eval=FALSE>>=
# Select ETF symbols for asset allocation
symbolv <- c("SPY", "VTI", "QQQ", "VEU", "EEM", "XLY", "XLP", 
"XLE", "XLF", "XLV", "XLI", "XLB", "XLK", "XLU", "VYM", "IVW", 
"IWB", "IWD", "IWF", "IEF", "TLT", "VNQ", "DBC", "GLD", "USO", 
"VXX", "SVXY", "MTUM", "IVE", "VLUE", "QUAL", "VTV", "USMV", "AIEQ")
# Setup code
etfenv <- new.env()  # New environment for data
# Boolean vector of symbols already downloaded
isdown <- symbolv %in% ls(etfenv)
      @
    \column{0.6\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Download data from Polygon using while loop
while (sum(!isdown) > 0) {
  for (symboln in symbolv[!isdown]) {
    cat("Processing:", symboln, "\n")
    tryCatch({  # With error handler
      # Download OHLC bars from Polygon into JSON format file
      urll <- paste0("https://api.polygon.io/v2/aggs/ticker/", symboln, "/range/1/", tspan, "/", startd, "/", todayd, "?adjusted=true&sort=asc&limit=50000&apiKey=", apikey)
      ohlc <- jsonlite::read_json(urll)
      # Extract list of prices from json object
      ohlc <- ohlc$results
      # Coerce from list to matrix
      ohlc <- lapply(ohlc, unlist)
      ohlc <- do.call(rbind, ohlc)
      # Coerce time from milliseconds to dates
      datev <- ohlc[, "t"]/1e3
      datev <- as.POSIXct(datev, origin="1970-01-01")
      datev <- as.Date(datev)
      # Coerce from matrix to xts
      ohlc <- ohlc[, c("o","h","l","c","v","vw")]
      colnames(ohlc) <- paste0(symboln, ".", c("Open", "High", "Low", "Close", "Volume", "VWAP"))
      ohlc <- xts::xts(ohlc, order.by=datev)
      # Save to environment
      assign(symboln, ohlc, envir=etfenv)
      Sys.sleep(1)
      },
    error={function(msg) print(paste0("Error handler: ", msg))},
    finally=print(paste0("Symbol = ", symboln))
    )  # end tryCatch
  }  # end for
  # Update vector of symbols already downloaded
  isdown <- symbolv %in% ls(etfenv)
}  # end while
save(etfenv, file="/Users/jerzy/Develop/lecture_slides/data/etf_data.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Calculating the Stock Alphas, Betas, and Other Performance Statistics}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package 
      \href{https://cran.r-project.org/web/packages/PerformanceAnalytics/index.html}{\emph{PerformanceAnalytics}} 
      contains functions for calculating risk and performance statistics, such as the \emph{variance}, \emph{skewness}, \emph{kurtosis}, \emph{beta}, \emph{alpha}, etc.
      \vskip1ex
      The function \texttt{PerformanceAnalytics::table.CAPM()} calculates the \emph{beta} $\beta$ and \emph{alpha} $\alpha$ values, the \emph{Treynor} ratio, and other performance statistics.
      \vskip1ex
      The function \texttt{PerformanceAnalytics::table.Stats()} calculates a data frame of risk and return statistics of the return distributions.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
# Extract Close prices
prices <- eapply(etfenv, quantmod::Cl)
prices <- do.call(cbind, prices)
# Drop ".Close" from colnames
colnames(prices) <- do.call(rbind, strsplit(colnames(prices), split="[.]"))[, 1]
# Calculate the log returns
retp <- xts::diff.xts(log(prices))
# Copy prices and returns into etfenv
etfenv$prices <- prices
etfenv$retp <- retp
# Copy symbolv into etfenv
etfenv$symbolv <- symbolv
# Calculate the risk-return statistics
riskstats <- PerformanceAnalytics::table.Stats(retp)
# Transpose the data frame
riskstats <- as.data.frame(t(riskstats))
# Add Name column
riskstats$Name <- rownames(riskstats)
# Copy riskstats into etfenv
etfenv$riskstats <- riskstats
# Calculate the beta, alpha, Treynor ratio, and other performance statistics
capmstats <- PerformanceAnalytics::table.CAPM(Ra=retp[, symbolv], 
                                               Rb=retp[, "VTI"], scale=252)
colv <- strsplit(colnames(capmstats), split=" ")
colv <- do.call(cbind, colv)[1, ]
colnames(capmstats) <- colv
capmstats <- t(capmstats)
capmstats <- capmstats[, -1]
colv <- colnames(capmstats)
whichv <- match(c("Annualized Alpha", "Information Ratio", "Treynor Ratio"), colv)
colv[whichv] <- c("Alpha", "Information", "Treynor")
colnames(capmstats) <- colv
capmstats <- capmstats[order(capmstats[, "Alpha"], decreasing=TRUE), ]
# Copy capmstats into etfenv
etfenv$capmstats <- capmstats
save(etfenv, file="/Users/jerzy/Develop/lecture_slides/data/etf_data.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Scraping \protect\emph{S\&P500} Stock Index Constituents From Websites}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{S\&P500} index constituents change over time, and \emph{Standard \& Poor's} replaces companies that have decreased in capitalization with ones that have increased.
      \vskip1ex
      The \emph{S\&P500} index may contain more than 500 stocks because some companies have several share classes of stock.
      \vskip1ex
      The \emph{S\&P500} index constituents may be scraped from websites like \href{https://en.wikipedia.org/wiki/List_of_S%26P_500_companies}{Wikipedia}, using dedicated packages.
      \vskip1ex
      The function \texttt{getURL()} from package \emph{RCurl} downloads the \emph{html} text data from an internet website \texttt{URL}.
      \vskip1ex
      The function \texttt{readHTMLTable()} from package \emph{XML} extracts tables from \emph{html} text data or from a remote \texttt{URL}, and returns them as a list of \emph{data frames} or matrices.
      \vskip1ex
      \texttt{readHTMLTable()} can't parse secure \texttt{URLs}, so they must first be downloaded using function \texttt{getURL()}, and then parsed using \texttt{readHTMLTable()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
library(RCurl)  # Load package RCurl
library(XML)  # Load package XML
# Download text data from URL
sp500 <- getURL(
  "https://en.wikipedia.org/wiki/List_of_S%26P500_companies")
# Extract tables from the text data
sp500 <- readHTMLTable(sp500)
str(sp500)
# Extract colnames of data frames
lapply(sp500, colnames)
# Extract S&P500 constituents
sp500 <- sp500[[1]]
head(sp500)
# Create valid R names from symbols containing "-" or "."characters
sp500$namev <- gsub("-", "_", sp500$Ticker)
sp500$namev <- gsub("[.]", "_", sp500$names)
# Write data frame of S&P500 constituents to CSV file
write.csv(sp500,
  file="/Users/jerzy/Develop/lecture_slides/data/sp500_Yahoo.csv",
  row.names=FALSE)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{S\&P500} Time Series Data From \protect\emph{Yahoo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Before time series data for the \emph{S\&P500} index constituents can be downloaded from \emph{Yahoo}, it's necessary to create valid names corresponding to symbols containing special characters like \texttt{"-"}.
      \vskip1ex
      The function \texttt{setSymbolLookup()} creates a lookup table for \emph{Yahoo} symbols, using valid names in \texttt{R}.
      \vskip1ex
      For example \emph{Yahoo} uses the symbol \texttt{"BRK-B"}, which isn't a valid name in \texttt{R}, but can be mapped to \texttt{"BRK\_B"}, using the function \texttt{setSymbolLookup()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
library(rutils)  # Load package rutils
# Load data frame of S&P500 constituents from CSV file
sp500 <- read.csv(file="/Users/jerzy/Develop/lecture_slides/data/sp500_Yahoo.csv")
# Register symbols corresponding to R names
for (indeks in 1:NROW(sp500)) {
  cat("processing: ", sp500$Ticker[indeks], "\n")
  setSymbolLookup(structure(
    list(list(name=sp500$Ticker[indeks])),
    names=sp500$names[indeks]))
}  # end for
sp500env <- new.env()  # new environment for data
# Remove all files (if necessary)
rm(list=ls(sp500env), envir=sp500env)
# Download data and copy it into environment
rutils::get_data(sp500$names,
   env_out=sp500env, startd="1990-01-01")
# Or download in loop
for (symboln in sp500$names) {
  cat("processing: ", symboln, "\n")
  rutils::get_data(symboln,
   env_out=sp500env, startd="1990-01-01")
}  # end for
save(sp500env, file="/Users/jerzy/Develop/lecture_slides/data/sp500.RData")
chart_Series(x=sp500env$BRKB["2016/"],
             TA="add_Vo()", name="BRK-B stock")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{FRED} Time Series Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{FRED} is a database of economic time series maintained by the Federal Reserve Bank of St. Louis:\\
      \hskip1em\url{http://research.stlouisfed.org/fred2/}
      \vskip1ex
      The function \texttt{getSymbols()} downloads time series data into the specified \emph{environment}.
      \vskip1ex
      \texttt{getSymbols()} can download \emph{FRED} data with the argument \texttt{"src"} set to \texttt{FRED}.
      \vskip1ex
      If the argument \texttt{"auto.assign"} is set to \texttt{FALSE}, then \texttt{getSymbols()} returns the data, instead of assigning it silently.
    \column{0.5\textwidth}
    \vspace{-1em}
      \includegraphics[width=0.5\paperwidth]{figure/fred_unemp_rate.png}
    \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Download U.S. unemployment rate data
unrate <- quantmod::getSymbols("UNRATE", 
   auto.assign=FALSE, src="FRED")
# Plot U.S. unemployment rate data
dygraphs::dygraph(unrate["1990/"], main="U.S. Unemployment Rate") %>%
  dyOptions(colors="blue", strokeWidth=2)
# Or
quantmod::chart_Series(unrate["1990/"], name="U.S. Unemployment Rate")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{The \protect\emph{Quandl} Database}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \href{https://www.quandl.com/}{Quandl} is a distributor of third party data, and offers several million financial, economic, and social datasets.
      \vskip1ex
      Much of the \href{https://www.quandl.com/}{Quandl} data is free, while premium data can be obtained under a temporary license.
      \vskip1ex
      \href{https://www.quandl.com/}{Quandl} provides online help and a guide to its datasets:\\
      \hskip1em\url{https://www.quandl.com/help/r}\\
      \hskip1em\url{https://www.quandl.com/browse}\\
      \hskip1em\url{https://www.quandl.com/blog/getting-started-with-the-quandl-api}\\
      \hskip1em\url{https://www.quandl.com/blog/stock-market-data-guide}
      \vskip1ex
      \href{https://www.quandl.com/}{Quandl} provides stock prices, stock fundamentals, financial ratios, zoo::indexes, options and volatility, earnings estimates, analyst ratings, etc.:\\
      \hskip1em\url{https://www.quandl.com/blog/api-for-stock-data}
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
install.packages("devtools")
library(devtools)
# Install package Quandl from github
install_github("quandl/R-package")
library(Quandl)  # Load package Quandl
# Register Quandl API key
Quandl.api_key("pVJi9Nv3V8CD3Js5s7Qx")
# Get short description
packageDescription("Quandl")
# Load help page
help(package="Quandl")
# Remove Quandl from search path
detach("package:Quandl")
      @
      \href{https://www.quandl.com/}{Quandl} has developed an \texttt{R} package called \emph{Quandl} that allows downloading data from \href{https://www.quandl.com/}{Quandl} directly into \texttt{R}.
      \vskip1ex
      To make more than 50 downloads a day, you need to register your \emph{Quandl} \emph{API key} using the function \texttt{Quandl.api\_key()},
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Time Series Data from \protect\emph{Quandl}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \href{https://www.quandl.com/}{Quandl} data can be downloaded directly into \texttt{R} using the function \texttt{Quandl()}.
      \vskip1ex
      The dots \texttt{"..."} argument of the \texttt{Quandl()} function accepts additional parameters to the \emph{Quandl API},\\
      \vskip1ex
      \href{https://www.quandl.com/}{Quandl} datasets have a unique \emph{Quandl code} in the format \texttt{"database/ticker"}, which can be found on the \href{https://www.quandl.com/}{Quandl} website for that dataset:\\
      \hskip1em\url{https://www.quandl.com/data/WIKI?keyword=aapl}
      \vskip1ex
      \emph{WIKI} is a user maintained free database of daily prices for 3,000 U.S. stocks,\\
      \hskip1em\url{https://www.quandl.com/data/WIKI}
      \vskip1ex
      \emph{SEC} is a free database of stock fundamentals extracted from \emph{SEC} \emph{10Q} and \emph{10K} filings (but not harmonized),\\
      \hskip1em\url{https://www.quandl.com/data/SEC}
      \vskip1ex
      \emph{RAYMOND} is a free database of harmonized stock fundamentals, based on the \emph{SEC} database,
      \hskip1em\url{https://www.quandl.com/data/RAYMOND-Raymond}
      \hskip1em\url{https://www.quandl.com/data/RAYMOND-Raymond?keyword=aapl}
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
library(rutils)  # Load package rutils
# Download EOD AAPL prices from WIKI free database
pricev <- Quandl(code="WIKI/AAPL",
  type="xts", startd="1990-01-01")
x11(width=14, height=7)
chart_Series(pricev["2016", 1:4], name="AAPL OHLC prices")
# Add trade volume in extra panel
add_TA(pricev["2016", 5])
# Download euro currency rates
pricev <- Quandl(code="BNP/USDEUR",
    startd="2013-01-01",
    endd="2013-12-01", type="xts")
# Download multiple time series
pricev <- Quandl(code=c("NSE/OIL", "WIKI/AAPL"),
    startd="2013-01-01", type="xts")
# Download AAPL gross profits
prof_it <- Quandl("RAYMOND/AAPL_GROSS_PROFIT_Q", type="xts")
chart_Series(prof_it, name="AAPL gross profits")
# Download Hurst time series
pricev <- Quandl(code="PE/AAPL_HURST",
    startd="2013-01-01", type="xts")
chart_Series(pricev["2016/", 1], name="AAPL Hurst")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Stock Index and Instrument Metadata on \protect\emph{Quandl}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Instrument metadata specifies properties of instruments, like its currency, contract size, tick value, delievery months, start date, etc.
      \vskip1ex
      \href{https://www.quandl.com/}{Quandl} provides instrument metadata for stock indices, futures, and currencies:\\
      \hskip1em\url{https://www.quandl.com/blog/useful-listv}
      \vskip1ex
      \href{https://www.quandl.com/}{Quandl} also provides constituents for stock indices, for example the \emph{S\&P500}, \emph{Dow Jones Industrial Average}, \emph{NASDAQ Composite}, \emph{FTSE 100}, etc.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
# Load S&P500 stock Quandl codes
sp500 <- read.csv(
  file="/Users/jerzy/Develop/lecture_slides/data/sp500_quandl.csv")
# Replace "-" with "_" in symbols
sp500$free_code <- gsub("-", "_", sp500$free_code)
head(sp500)
# vector of symbols in sp500 frame
tickers <- gsub("-", "_", sp500$ticker)
# Or
tickers <- matrix(unlist(
  strsplit(sp500$free_code, split="/"),
  use.names=FALSE), ncol=2, byrow=TRUE)[, 2]
# Or
tickers <- do_call_rbind(
  strsplit(sp500$free_code, split="/"))[, 2]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Multiple Time Series from \protect\emph{Quandl}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Time series data for a portfolio of stocks can be downloaded by performing a loop over the function \texttt{Quandl()} from package \href{https://cran.r-project.org/web/packages/Quandl/index.html}{Quandl}.
      \vskip1ex
      The \texttt{assign()} function assigns a value to an object in a specified \emph{environment}, by referencing it using a character string (name).
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
sp500env <- new.env()  # new environment for data
# Remove all files (if necessary)
rm(list=ls(sp500env), envir=sp500env)
# Boolean vector of symbols already downloaded
isdown <- tickers %in% ls(sp500env)
# Download data and copy it into environment
for (ticker in tickers[!isdown]) {
  cat("processing: ", ticker, "\n")
  datav <- Quandl(code=paste0("WIKI/", ticker),
                  startd="1990-01-01", type="xts")[, -(1:7)]
  colnames(datav) <- paste(ticker,
    c("Open", "High", "Low", "Close", "Volume"), sep=".")
  assign(ticker, datav, envir=sp500env)
}  # end for
save(sp500env, file="/Users/jerzy/Develop/lecture_slides/data/sp500.RData")
chart_Series(x=sp500env$XOM["2016/"], TA="add_Vo()", name="XOM stock")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading Futures Time Series from \protect\emph{Quandl}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \href{https://www.quandl.com/}{Quandl} provides the
      \href{https://www.quandl.com/data/CHRIS-Wiki-Continuous-Futures}{Wiki CHRIS Database}
      of time series of prices for \texttt{600} different futures contracts.
      \vskip1ex
      The
      \href{https://www.quandl.com/data/CHRIS-Wiki-Continuous-Futures}{Wiki CHRIS Database} contains daily \emph{OHLC} prices for continuous futures contracts.
      \vskip1ex
      A continuous futures contract is a time series of prices obtained by chaining together prices from consecutive futures contracts.
      \vskip1ex
      The data is curated by the \href{https://www.quandl.com/}{Quandl} community from data provided by the \emph{CME}, \emph{ICE}, \emph{LIFFE}, and other exchanges.
      \vskip1ex
      The \emph{Quandl codes} are specified as \texttt{CHRIS/\{EXCHANGE\}\_\{CODE\}\{DEPTH\}}, where \texttt{\{DEPTH\}} is the depth of the chained contract.
      \vskip1ex
      The chained front month contracts have depth $1$, the back month contracts have depth $2$, etc.
      \vskip1ex
      The continuous front and back month contracts allow building continuous futures curves.
      \vskip1ex
      \href{https://www.quandl.com/}{Quandl} data can be downloaded directly into \texttt{R} using the function \texttt{Quandl()}.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)
library(Quandl)
# Register Quandl API key
Quandl.api_key("pVJi9Nv3V8CD3Js5s7Qx")
# Download E-mini S&P500 futures prices
pricev <- Quandl(code="CHRIS/CME_ES1",
  type="xts", startd="1990-01-01")
pricev <- pricev[, c("Open", "High", "Low", "Last", "Volume")]
colnames(pricev)[4] <- "Close"
# Plot the prices
x11(width=5, height=4)  # Open x11 for plotting
chart_Series(x=pricev["2008-06/2009-06"],
             TA="add_Vo()", name="S&P500 Futures")
# Plot dygraph
dygraphs::dygraph(pricev["2008-06/2009-06", -5],
  main="S&P500 Futures") %>%
  dyCandlestick()
      @
      \vspace{-1em}
      For example, the \emph{Quandl code} for the continuous \emph{E-mini S\&P500} front month futures is \texttt{CHRIS/CME\_ES1}, while for the back month it's \texttt{CHRIS/CME\_ES2}, for the second back month it's \texttt{CHRIS/CME\_ES3}, etc.
      \vskip1ex
      The \emph{Quandl code} for the \emph{E-mini Oil} futures is \texttt{CHRIS/CME\_QM1}, for the \emph{E-mini euro FX} futures is \texttt{CHRIS/CME\_E71}, etc.
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{VIX} Futures Files from CBOE}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The CFE (CBOE Futures Exchange) provides daily \href{https://markets.cboe.com/us/futures/market_statistics/historical_data/}{CBOE Historical Data for Volatility Futures}, including the \emph{VIX} futures.
      \vskip1ex
      The CBOE data incudes \emph{OHLC} prices and also the \emph{settlement} price (in column \texttt{"Settle"}).
      \vskip1ex
      The \emph{settlement} price is usually defined as the weighted average price (\emph{WAP}) or the midpoint price, and is different from the \emph{Close} price.
      \vskip1ex
      The \emph{settlement} price is used for calculating the daily \emph{mark to market} (value) of the futures contract.
      \vskip1ex
      Futures exchanges require that counterparties exchange (settle) the \emph{mark to market} value of the futures contract daily, to reduce counterparty default risk.
      \vskip1ex
      The function \texttt{download.file()} downloads files from the internet.
      \vskip1ex
      The function \texttt{tryCatch()} executes functions and expressions, and handles any \emph{exception conditions} produced when they are evaluated.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Read CBOE futures expiration dates
datev <- read.csv(file="/Users/jerzy/Develop/lecture_slides/data/futures_expiration_dates_codes.csv",
  row.names=1)
dirn <- "/Users/jerzy/Develop/data/vix_data"
dir.create(dirn)
symbolv <- rownames(datev)
filens <- file.path(dirn, paste0(symbolv, ".csv"))
log_file <- file.path(dirn, "log_file.txt")
cboe_url <- "https://markets.cboe.com/us/futures/market_statistics/historical_data/products/csv/VX/"
urls <- paste0(cboe_url, datev[, 1])
# Download files in loop
for (it in seq_along(urls)) {
    tryCatch(  # Warning and error handler
        download.file(urls[it],
                destfile=filens[it], quiet=TRUE),
      # Warning handler captures warning condition
      warning=function(msg) {
        cat(paste0("Warning handler: ", msg, "\n"), file=log_file, append=TRUE)
      },  # end warning handler
      # Error handler captures error condition
      error=function(msg) {
        cat(paste0("Error handler: ", msg, "\n"), append=TRUE)
      },  # end error handler
      finally=cat(paste0("Processing file name = ", filens[it], "\n"), append=TRUE)
    )  # end tryCatch
}  # end for
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{VIX} Futures Data Into an Environment}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{quantmod::getSymbols()} with the parameter \texttt{src="cfe"} downloads CFE data into the specified \emph{environment}. (But this requires first loading the package \emph{qmao}.)
      \vskip1ex
      Currently \texttt{quantmod::getSymbols()} doesn't download the most recent data.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Create new environment for data
vixenv <- new.env()
# Download VIX data for the months 6, 7, and 8 in 2018
library(qmao)
quantmod::getSymbols("VX", Months=1:12,
  Years=2018, src="cfe", auto.assign=TRUE, env=vixenv)
# Or
qmao::getSymbols.cfe(Symbols="VX",
  Months=6:8, Years=2018, env=vixenv,
  verbose=FALSE, auto.assign=TRUE)
# Calculate the classes of all the objects
# In the environment vixenv
unlist(eapply(vixenv, function(x) {class(x)[1]}))
class(vixenv$VX_M18)
colnames(vixenv$VX_M18)
# Save the data to a binary file called "vix_cboe.RData".
save(vixenv,
  file="/Users/jerzy/Develop/data/vix_data/vix_cboe.RData")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\section{High Frequency and Intraday Stock Prices}


%%%%%%%%%%%%%%%
\subsection{Trade and Quote (\protect\emph{TAQ}) Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{High frequency data} is typically formatted as either Trade and Quote (\emph{TAQ}) data, or \emph{Open-High-Low-Close} (\emph{OHLC}) data.
      \vskip1ex
      Trade and Quote (\emph{TAQ}) data contains intraday \emph{trades} and \emph{quotes} on exchange-traded stocks and futures.
      \vskip1ex
      \emph{TAQ} data is often called \emph{tick data}, with a \emph{tick} being a row of data containing new \emph{trades} or \emph{quotes}.
      \vskip1ex
      The \emph{TAQ} data is spaced irregularly in time, with data recorded each time a new trade or quote arrives.
      \vskip1ex
      Each row of \emph{TAQ} data may contain the quote and trade prices, and the corresponding quote size or trade volume:
      \emph{Bid.Price, Bid.Size, Ask.Price, Ask.Size, Trade.Price, Volume}.
      \vskip1ex
      \emph{TAQ} data is often split into \emph{trade} data and \emph{quote} data.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Load package HighFreq
library(HighFreq)
# Or load the high frequency data file directly:
# symbolv <- load("/Users/jerzy/Develop/R/HighFreq/data/hf_data.RData")
head(HighFreq::SPY_TAQ)
head(HighFreq::SPY)
tail(HighFreq::SPY)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Downloading \protect\emph{TAQ} Data From \protect\emph{WRDS}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{TAQ} data can be downloaded from the
      \href{https://wrds-web.wharton.upenn.edu/wrds/ds/taq/ctm/index.cfm}{\emph{WRDS TAQ}} web page.
      \vskip1ex
      The \emph{TAQ} data are at millisecond frequency, and are \emph{consolidated} (combined) from the New York Stock Exchange \emph{NYSE} and other exchanges.
      \vskip1ex
      The
      \href{https://wrds-web.wharton.upenn.edu/wrds/ds/taq/ctm/index.cfm}{\emph{WRDS TAQ}}
      web page provides separately \emph{trades} data and separately \emph{quotes} data.
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/wrds_taq_data.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Reading \protect\emph{TAQ} Data From \texttt{.csv} Files}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Trade and Quote (\emph{TAQ}) data stored in \texttt{.csv} files can be very large, so it's better to read it using the function \texttt{data.table::fread()} which is much faster than the function \texttt{read.csv()}.
      \vskip1ex
      Each \emph{trade} or \emph{quote} contributes a \emph{tick} (row) of data, and the number of ticks can be very large (hundred of thousands per day, or more).
      \vskip1ex
      The function \texttt{strptime()} coerces \texttt{character} strings representing the date and time into \texttt{POSIXlt} \emph{date-time} objects.
      \vskip1ex
      The argument \texttt{format="\%H:\%M:\%OS"} allows the parsing of fractional seconds, for example \texttt{"15:59:59.989847074"}.
      \vskip1ex
      The function \texttt{as.POSIXct()} coerces objects into \texttt{POSIXct} \emph{date-time} objects, with a \texttt{numeric} value representing the \emph{moment of time} in seconds.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
library(HighFreq)
# Read TAQ trade data from csv file
taq <- data.table::fread(file="/Users/jerzy/Develop/lecture_slides/data/xlk_tick_trades_20200316.csv")
# Inspect the TAQ data in data.table format
taq
class(taq)
colnames(taq)
sapply(taq, class)
symboln <- taq$SYM_ROOT[1]
# Create date-time index
datev <- paste(taq$DATE, taq$TIME_M)
# Coerce date-time index to POSIXlt
datev <- strptime(datev, "%Y%m%d %H:%M:%OS")
class(datev)
# Display more significant digits
# options("digits")
options(digits=20, digits.secs=10)
last(datev)
unclass(last(datev))
as.numeric(last(datev))
# Coerce date-time index to POSIXct
datev <- as.POSIXct(datev)
class(datev)
last(datev)
unclass(last(datev))
as.numeric(last(datev))
# Calculate the number of seconds
as.numeric(last(datev)) - as.numeric(first(datev))
# Calculate the number of ticks per second
NROW(taq)/(6.5*3600)
# Select TAQ data columns
taq <- taq[, .(price=PRICE, volume=SIZE)]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Trading Volumes in High Frequency Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The trading volumes represent the number of shares traded at a given price.  
      \vskip1ex
      The histogram of the trading volumes shows that the highest frequencies of trades are for \texttt{100} shares and for round lots (trades that are multiples of \texttt{100} shares.)
      \vskip1ex
      There are also significant frequencies for \emph{odd lots}, with small volumes of less than \texttt{100} shares.
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.4\paperwidth]{figure/hf_ticks_volume_histogram.png}
      <<echo=TRUE,eval=FALSE>>=
# Coerce trade ticks to xts series
xlk <- xts::xts(taq[, .(price, volume)], datev)
colnames(xlk) <- c("price", "volume")
save(xlk, file="/Users/jerzy/Develop/data/xlk_tick_trades_20200316.RData")
# Plot histogram of the trading volumes
hist(xlk$volume, main="Histogram of XLK Trading Volumes", 
     breaks=1e5, xlim=c(1, 400), xlab="number of shares")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Microstructure Noise in High Frequency Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      High frequency data contains \emph{microstructure noise} in the form of \emph{price spikes} and the \emph{bid-ask bounce}.
      \vskip1ex
      \emph{Price spikes} are single ticks with prices far away from the average.
      \vskip1ex
      \emph{Price spikes} are often caused by data collection errors, but sometimes they represent actual trades with very large lot (trade) sizes.
      \vskip1ex
      The \emph{bid-ask bounce} is the bouncing of traded prices between the bid and ask prices.
      \vskip1ex
      The \emph{bid-ask bounce} creates an illusion of rapidly changing prices, while in reality the mid price is unchanged.
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_raw.png}
      % \vspace{-1em}
      % \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_bounce.png}
      <<echo=TRUE,eval=FALSE>>=
# Plot dygraph
dygraphs::dygraph(xlk$price, main="XLK Intraday Prices for 2020-03-16") %>%
  dyOptions(colors="blue", strokeWidth=1)
# Plot in x11 window
x11(width=6, height=5)
quantmod::chart_Series(x=xlk$price, name="XLK Intraday Prices for 2020-03-16")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{The Bid-ask Bounce of High Frequency Prices}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{bid-ask bounce} is the bouncing of traded prices between the bid and ask prices.
      \vskip1ex
      The \emph{bid-ask bounce} is prominent at very high frequency time scales or in periods of low volatility.
      \vskip1ex
      The \emph{bid-ask bounce} creates an illusion of rapidly changing prices, while in fact the mid price is constant.
      \vskip1ex
      The \emph{bid-ask bounce} inflates the estimates of realized volatility, above the actual volatility.
      \vskip1ex
      The \emph{bid-ask bounce} creates the appearance of mean reversion (negative autocorrelation), that isn't tradeable for most traders.
      <<echo=TRUE,eval=FALSE>>=
pricev <- read.zoo(file="/Users/jerzy/Develop/lecture_slides/data/bid_ask_bounce.csv",
  header=TRUE, sep=",")
pricev <- as.xts(pricev)
dygraphs::dygraph(pricev$Close, 
  main="S&P500 Futures Prices Bid-Ask Bounce") %>%
  dyOptions(colors="blue", strokeWidth=2)
      @
    \column{0.5\textwidth}
      \includegraphics[width=0.45\paperwidth]{figure/hf_futures_bounce.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Price Spikes And Trading Volumes in High Frequency Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The number of the \emph{price spikes} depends on the level of trading volumes, with the number decreasing with higher trading volumes.
      \vskip1ex
      The number of price spikes is lower for trade prices with larger trading volumes.
    \column{0.5\textwidth}
      \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_biglots.png}
      % \vspace{-1em}
      % \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_bounce.png}
      <<echo=TRUE,eval=FALSE>>=
# Plot dygraph of trade prices of at least 100 shares
dygraphs::dygraph(xlk$price[xlk$volume >= 100, ],
  main="XLK Prices for Trades of At Least 100 Shares") %>%
  dyOptions(colors="blue", strokeWidth=1)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Removing Odd Lot Trades From \protect\emph{TAQ} Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The trading volumes represent the number of shares traded at a given price.  
      \vskip1ex
      The histogram of the trading volumes shows that the highest frequencies are for \texttt{100} shares and for round lots (trades that are multiples of \texttt{100} shares.)
      \vskip1ex
      There are also significant frequencies for \emph{odd lots}, with small volumes of less than \texttt{100} shares.
      \vskip1ex
      The \emph{odd lot} ticks are often removed to reduce the size of the \emph{TAQ} data.
      \vskip1ex
      Selecting only the large lot trades reduces microstructure noise (price spikes, bid-ask bounce) in high frequency data.
      <<echo=TRUE,eval=FALSE>>=
# Select the large trade lots of at least 100 shares
dim(taq)
tickb <- taq[taq$volume >= 100]
dim(tickb)
# Number of large lot ticks per second
NROW(tickb)/(6.5*3600)
# Plot histogram of the trading volumes
hist(tickb$volume, main="Histogram of XLK Trading Volumes", 
     breaks=100000, xlim=c(1, 400), xlab="number of shares")
# Save trade ticks with large lots
data.table::fwrite(tickb, file="/Users/jerzy/Develop/data/xlk_tick_trades_20200316_biglots.csv")
# Coerce trade prices to xts
xlkb <- xts::xts(tickb[, .(price, volume)], tickb$index)
colnames(xlkb) <- c("price", "volume")
      @
    \column{0.5\textwidth}
      \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_biglots.png}
      <<echo=TRUE,eval=FALSE>>=
# Plot dygraph of the large lots
dygraphs::dygraph(xlkb$price, 
  main="XLK Prices for Trades of At Least 100 Shares") %>%
  dyOptions(colors="blue", strokeWidth=1)
# Plot the large lots
x11(width=6, height=5)
quantmod::chart_Series(x=xlk$price,
  name="XLK Trade Ticks for 2020-03-16 (large lots only)")
      @
  \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Scrubbing Bad Market Data}


%%%%%%%%%%%%%%%
\subsection{The Hampel Filter For Filtering Price Spikes}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Price spikes in high frequency data can be identified using a \emph{Hampel filter}.
      \vskip1ex
      The \emph{z-scores} are equal to the prices minus the median of the prices, divided by the median absolute deviation (\emph{MAD}) of prices:
      \begin{displaymath}
        z_i = \frac{p_i - \operatorname{median}(\mathbf{p})}{\operatorname{MAD}}
      \end{displaymath}
      If the absolute value of the \emph{z-score} exceeds the \emph{threshold value} then it's classified as \emph{bad data}, and it can be removed or replaced.
      <<echo=TRUE,eval=FALSE>>=
# Calculate the centered Hampel filter to remove bad prices
lookb <- 71 # Look-back interval
halfb <- lookb %/% 2 # Half-back interval
pricev <- xlk$price
# Calculate the trailing median and MAD
medianv <- HighFreq::roll_mean(pricev, lookb=lookb, method="nonparametric")
colnames(medianv) <- c("median")
madv <- HighFreq::roll_var(pricev, lookb=lookb, method="nonparametric")
# madv <- TTR::runMAD(pricev, n=lookb)
# Center the median and the MAD
medianv <- rutils::lagit(medianv, lagg=(-halfb), pad_zeros=FALSE)
madv <- rutils::lagit(madv, lagg=(-halfb), pad_zeros=FALSE)
# Calculate the Z-scores
zscores <- ifelse(madv > 0, (pricev - medianv)/madv, 0)
# Z-scores have very fat tails
range(zscores); mad(zscores)
madz <- mad(zscores[abs(zscores) > 0])
hist(zscores, breaks=50000, xlim=c(-2*madz, 2*madz))
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_scrubbed.png}
      <<echo=TRUE,eval=FALSE>>=
# Define discrimination threshold value
threshv <- 6*madz
# Identify good prices with small z-scores
isgood <- (abs(zscores) < threshv)
# Calculate the number of bad prices
sum(!isgood)
# Overwrite bad prices and calculate time series of scrubbed prices
priceg <- pricev
priceg[!isgood] <- NA
priceg <- zoo::na.locf(priceg)
# Plot dygraph of the scrubbed prices
dygraphs::dygraph(priceg, main="Scrubbed XLK Intraday Prices") %>%
  dyOptions(colors="blue", strokeWidth=1)
# Plot using chart_Series()
x11(width=6, height=5)
quantmod::chart_Series(x=priceg,
  name="Clean XLK Intraday Prices for 2020-03-16")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Classifying Data Outliers Using the Hampel Filter}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The Hampel filter is a \emph{classifier} which classifies the prices as either good or bad data points.
      \vskip1ex
      In order to measure the performance of the Hampel filter, we add price spikes to the clean prices, to see how accurately they're classified. 
      \vskip1ex
      Let the \emph{null hypothesis} be that the given price is a good data point.
      \vskip1ex
      A positive result corresponds to rejecting the \emph{null hypothesis}, while a negative result corresponds to accepting the \emph{null hypothesis}.
      \vskip1ex
      The classifications are subject to two different types of errors: \emph{type I} and \emph{type II} errors.
      \vskip1ex
      A \emph{type I} error is the incorrect rejection of a \texttt{TRUE} \emph{null hypothesis} (i.e. a "false positive"), when good data is classified as bad.
      \vskip1ex
      A \emph{type II} error is the incorrect acceptance of a \texttt{FALSE} \emph{null hypothesis} (i.e. a "false negative"), when bad data is classified as good.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Add 200 random price spikes to the clean prices
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
nspikes <- 200
nrows <- NROW(priceg)
ispike <- logical(nrows)
ispike[sample(x=nrows, size=nspikes)] <- TRUE
priceb <- priceg
priceb[ispike] <- priceb[ispike]*
  sample(c(0.999, 1.001), size=nspikes, replace=TRUE)
# Plot the bad prices and their medians
medianv <- HighFreq::roll_mean(priceb, lookb=lookb, method="nonparametric")
pricem <- cbind(priceb, medianv)
colnames(pricem) <- c("prices with spikes", "median")
dygraphs::dygraph(pricem, main="XLK Prices With Spikes") %>%
  dyOptions(colors=c("red", "blue"))
# Calculate the z-scores
madv <- HighFreq::roll_var(priceb, lookb=lookb, method="nonparametric")
zscores <- ifelse(madv > 0, (priceb - medianv)/madv, 0)
# Z-scores have very fat tails
range(zscores); mad(zscores)
madz <- mad(zscores[abs(zscores) > 0])
hist(zscores, breaks=10000, xlim=c(-4*madz, 4*madz))
# Identify good prices with small z-scores
threshv <- 3*madz
isgood <- (abs(zscores) < threshv)
# Calculate the number of bad prices
sum(!isgood)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Confusion Matrix of a Binary Classification Model}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      A \emph{binary classification model} categorizes cases based on its forecasts whether the \emph{null hypothesis} is \texttt{TRUE} or \texttt{FALSE}.
      \vskip1ex
      The confusion matrix summarizes the performance of a classification model on a set of test data for which the actual values of the \emph{null hypothesis} are known.
      \vskip1ex
      \newcommand\MyBox[2]{
        \fbox{\lower0.75cm
          \vbox to 1.2cm{\vfil
            \hbox to 1.7cm{\parbox{\textwidth}{#1\\#2}}
            \vfil}
        }
      }
      \renewcommand\arraystretch{0.3}
      \setlength\tabcolsep{0pt}
      {\tiny
      \begin{tabular}{c >{\bfseries}r @{\hspace{0.5em}}c @{\hspace{0.4em}}c @{\hspace{0.5em}}l}
      \multirow{10}{*}{\parbox{0.5cm}{\bfseries Actual}} &
      & \multicolumn{2}{c}{\bfseries Forecast} & \\
      & & \bfseries Null is FALSE & \bfseries Null is TRUE \\
      & Null is FALSE & \MyBox{True Positive}{(sensitivity)} & \MyBox{False Negative}{(type II error)} \\[2.4em]
      & Null is TRUE & \MyBox{False Positive}{(type I error)} & \MyBox{True Negative}{(specificity)}
      \end{tabular}}
      <<echo=TRUE,eval=FALSE>>=
# Calculate the confusion matrix
table(actual=!ispike, forecast=isgood)
sum(!isgood)
# FALSE positive (type I error)
sum(!ispike & !isgood)
# FALSE negative (type II error)
sum(ispike & isgood)
      @
    \column{0.5\textwidth}
      Let the \emph{null hypothesis} be that the given price is a good data point.
      \vskip1ex
      The \emph{true positive} rate (known as the \emph{sensitivity}) is the fraction of \texttt{FALSE} \emph{null hypothesis} cases that are correctly classified as \texttt{FALSE}.
      \vskip1ex
      The \emph{false negative} rate is the fraction of \texttt{FALSE} \emph{null hypothesis} cases that are incorrectly classified as \texttt{TRUE} (\emph{type II} error).
      \vskip1ex
      The sum of the \emph{true positive} plus the \emph{false negative} rate is equal to $1$.
      \vskip1ex
      The \emph{true negative} rate (known as the \emph{specificity}) is the fraction of \texttt{TRUE} \emph{null hypothesis} cases that are correctly classified as \texttt{TRUE}.
      \vskip1ex
      The \emph{false positive} rate is the fraction of \texttt{TRUE} \emph{null hypothesis} cases that are incorrectly classified as \texttt{FALSE} (\emph{type I} error).
      \vskip1ex
      The sum of the \emph{true negative} plus the \emph{false positive} rate is equal to $1$.
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Receiver Operating Characteristic (ROC) Curve}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{ROC curve} is the plot of the \emph{true positive} rate, as a function of the \emph{false positive} rate, and illustrates the performance of a binary classifier.
      \vskip1ex
      The area under the \emph{ROC curve} (AUC) measures the classification ability of a binary classifier.
      <<echo=TRUE,eval=FALSE>>=
# Confusion matrix as function of threshold
confun <- function(actualv, zscores, threshv) {
    confmat <- table(actualv, (abs(zscores) < threshv))
    confmat <- confmat / rowSums(confmat)
    c(typeI=confmat[2, 1], typeII=confmat[1, 2])
}  # end confun
confun(!ispike, zscores, threshv=threshv)
# Define vector of discrimination thresholds
threshv <- madz*seq(from=0.1, to=3.0, by=0.05)/2
# Calculate the error rates
errorr <- sapply(threshv, confun, actualv=!ispike, zscores=zscores)
errorr <- t(errorr)
rownames(errorr) <- threshv
errorr <- rbind(c(1, 0), errorr)
errorr <- rbind(errorr, c(0, 1))
# Calculate the area under the ROC curve (AUC)
truepos <- (1 - errorr[, "typeII"])
truepos <- (truepos + rutils::lagit(truepos))/2
falsepos <- rutils::diffit(errorr[, "typeI"])
abs(sum(truepos*falsepos))
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hampel_roc.png}
      <<echo=TRUE,eval=FALSE>>=
# Plot ROC curve for Hampel classifier
plot(x=errorr[, "typeI"], y=1-errorr[, "typeII"],
     xlab="FALSE positive rate", ylab="TRUE positive rate",
     xlim=c(0, 1), ylim=c(0, 1),
     main="ROC Curve for Hampel Classifier",
     type="l", lwd=3, col="blue")
abline(a=0.0, b=1.0, lwd=3, col="orange")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Filtering Bad Data From Daily Stock Prices}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Daily stock prices can also contain bad data points consisting of mostly single, isolated spikes in prices.
      \vskip1ex
      The number of false positives may be too high, so the Hampel filter parameters (the look-back interval and the threshold) need adjustment.
      \vskip1ex
      For example, the \emph{VXX} has only one bad price (on \texttt{2010-11-08}), but the Hampel filter identifies many more than that (which are false positives).
      <<echo=TRUE,eval=FALSE>>=
# Load log VXX prices
load("/Users/jerzy/Develop/lecture_slides/data/pricevxx.RData")
nrows <- NROW(pricev)
# Calculate the centered Hampel filter for VXX
lookb <- 7 # Look-back interval
halfb <- lookb %/% 2 # Half-back interval
medianv <- HighFreq::roll_mean(pricev, lookb=lookb, method="nonparametric")
medianv <- rutils::lagit(medianv, lagg=(-halfb), pad_zeros=FALSE)
madv <- HighFreq::roll_var(pricev, lookb=lookb, method="nonparametric")
madv <- rutils::lagit(madv, lagg=(-halfb), pad_zeros=FALSE)
zscores <- ifelse(madv > 0, (pricev - medianv)/madv, 0)
range(zscores); mad(zscores)
madz <- mad(zscores[abs(zscores) > 0])
hist(zscores, breaks=100, xlim=c(-3*madz, 3*madz))
# Define discrimination threshold value
threshv <- 9*madz
# Calculate the good prices
isgood <- (abs(zscores) < threshv)
sum(!isgood)
# Dates of the bad prices
zoo::index(pricev[!isgood])
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hampel_stock_falsepos.png}
      <<echo=TRUE,eval=FALSE>>=
# Calculate the false positives
falsep <- !isgood
falsep[which(zoo::index(pricev) == as.Date("2010-11-08"))] <- FALSE
# Plot dygraph of the prices with bad prices
datam <- cbind(pricev, zscores)
colnames(datam)[2] <- "ZScores"
colv <- colnames(datam)
dygraphs::dygraph(datam, main="VXX Prices With Z-Scores and False Positives") %>%
  dyAxis("y", label=colv[1], independentTicks=TRUE) %>%
  dyAxis("y2", label=colv[2], independentTicks=TRUE) %>%
  dySeries(name=colv[1], axis="y", strokeWidth=1, col="blue") %>%
  dySeries(name=colv[2], axis="y2", strokeWidth=1, col="red") %>%
  dyEvent(zoo::index(pricev[falsep]), label=rep("false", sum(falsep)), strokePattern="solid", color="red") %>%
  dyEvent(zoo::index(pricev["2010-11-08"]), label="true", strokePattern="solid", color="green")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Scrubbing Bad Stock Prices}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Bad stock prices can be scrubbed (replaced) with the previous good price.  
      \vskip1ex
      But it's incorrect to replace bad prices with the average of the previous good price and the next good price, since that would cause data snooping.
      <<echo=TRUE,eval=FALSE>>=
# Replace bad stock prices with the previous good prices
priceg <- pricev
priceg[!isgood] <- NA
priceg <- zoo::na.locf(priceg)
# Calculate the Z-scores
medianv <- HighFreq::roll_mean(priceg, lookb=lookb, method="nonparametric")
medianv <- rutils::lagit(medianv, lagg=(-halfb), pad_zeros=FALSE)
madv <- HighFreq::roll_var(priceg, lookb=lookb, method="nonparametric")
madv <- rutils::lagit(madv, lagg=(-halfb), pad_zeros=FALSE)
zscores <- ifelse(madv > 0, (priceg - medianv)/madv, 0)
madz <- mad(zscores[abs(zscores) > 0])
# Calculate the number of bad prices
threshv <- 9*madz
isgood <- (abs(zscores) < threshv)
sum(!isgood)
zoo::index(priceg[!isgood])
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hampel_stock_scrubbed.png}
      <<echo=TRUE,eval=FALSE>>=
# Calculate the false positives
falsep <- !isgood
falsep[which(zoo::index(pricev) == as.Date("2010-11-08"))] <- FALSE
# Plot dygraph of the prices with bad prices
dygraphs::dygraph(priceg, main="Scrubbed VXX Prices With False Positives") %>%
  dyEvent(zoo::index(priceg[falsep]), label=rep("false", sum(falsep)), strokePattern="solid", color="red") %>%
  dyOptions(colors="blue", strokeWidth=1)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{ROC Curve for Daily Hampel Classifier}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      In order to measure the performance of the Hampel filter, we add price spikes to the clean prices, to see how accurately they're classified. 
      \vskip1ex
      The performance of the Hampel noise classification model depends on the length of the look-back time interval.
      \vskip1ex
      The optimal \emph{look-back interval} and \emph{threshold value} can be determined using \emph{cross-validation}. 
      <<echo=TRUE,eval=FALSE>>=
# Add 200 random price spikes to the clean prices
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
nspikes <- 200
ispike <- logical(nrows)
ispike[sample(x=nrows, size=nspikes)] <- TRUE
priceb <- priceg
priceb[ispike] <- priceb[ispike]*
  sample(c(0.99, 1.01), size=nspikes, replace=TRUE)
# Calculate the Z-scores
medianv <- HighFreq::roll_mean(priceb, lookb=lookb, method="nonparametric")
medianv <- rutils::lagit(medianv, lagg=(-halfb), pad_zeros=FALSE)
madv <- HighFreq::roll_var(priceb, lookb=lookb, method="nonparametric")
madv <- rutils::lagit(madv, lagg=(-halfb), pad_zeros=FALSE)
zscores <- ifelse(madv > 0, (priceb - medianv)/madv, 0)
madz <- mad(zscores[abs(zscores) > 0])
# Define vector of discrimination thresholds
threshv <- madz*seq(from=0.1, to=3.0, by=0.05)/2
# Calculate the error rates
errorr <- sapply(threshv, confun, actualv=!ispike, zscores=zscores)
errorr <- t(errorr)
rownames(errorr) <- threshv
errorr <- rbind(c(1, 0), errorr)
errorr <- rbind(errorr, c(0, 1))
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hampel_roc_daily.png}
      <<echo=TRUE,eval=FALSE>>=
# Calculate the area under the ROC curve (AUC)
truepos <- (1 - errorr[, "typeII"])
truepos <- (truepos + rutils::lagit(truepos))/2
falsepos <- rutils::diffit(errorr[, "typeI"])
abs(sum(truepos*falsepos))
# Plot ROC curve for Hampel classifier
plot(x=errorr[, "typeI"], y=1-errorr[, "typeII"],
     xlab="FALSE positive rate", ylab="TRUE positive rate",
     xlim=c(0, 1), ylim=c(0, 1),
     main="ROC Curve for Daily Hampel Classifier",
     type="l", lwd=3, col="blue")
abline(a=0.0, b=1.0, lwd=3, col="orange")
      @
  \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Modeling of Intraday Stock Prices}


%%%%%%%%%%%%%%%
\subsection{Aggregating \protect\emph{TAQ} Data to \protect\emph{OHLC}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{data table} columns can be \emph{aggregated} over categories (factors) defined by one or more columns passed to the \texttt{"by"} argument.
      \vskip1ex
      Multiple \emph{data table} columns can be referenced by passing a list of names specified by the dot \texttt{.()} operator.
      \vskip1ex
      The function \texttt{round.POSIXt()} rounds date-time objects to seconds, minutes, hours, days, months or years.
      \vskip1ex
      The function \texttt{as.POSIXct()} coerces objects to class \texttt{POSIXct}.
      <<echo=TRUE,eval=FALSE>>=
# Round time index to seconds
tickg[, zoo::index := as.POSIXct(round.POSIXt(index, "secs"))]
# Aggregate to OHLC by seconds
ohlc <- tickg[, .(open=first(price), high=max(price), low=min(price), close=last(price), volume=sum(volume)), by=index]
# Round time index to minutes
tickg[, zoo::index := as.POSIXct(round.POSIXt(index, "mins"))]
# Aggregate to OHLC by minutes
ohlc <- tickg[, .(open=first(price), high=max(price), low=min(price), close=last(price), volume=sum(volume)), by=index]
      @
    \column{0.5\textwidth}
      \includegraphics[width=0.45\paperwidth]{figure/hf_ticks_ohlc.png}
      <<echo=TRUE,eval=FALSE>>=
# Coerce OHLC prices to xts
ohlc <- xts::xts(ohlc[, -"index"], ohlc$index)
# Plot dygraph of the OHLC prices
dygraphs::dygraph(ohlc[, -5], main="XLK Trade Ticks for 2020-03-16 (OHLC)") %>%
  dyCandlestick()
# Plot the OHLC prices
x11(width=6, height=5)
quantmod::chart_Series(x=ohlc, TA="add_Vo()",
  name="XLK Trade Ticks for 2020-03-16 (OHLC)")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Open-High-Low-Close (\protect\emph{OHLC}) Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Open-High-Low-Close} (\emph{OHLC}) data contains intraday trade prices and trade volumes.
      \vskip1ex
      \emph{OHLC} data is evenly spaced in time, with each row containing the \emph{Open, High, Low, Close} prices, and the trade \emph{Volume}, recorded over the past time interval (called a \emph{bar} of data).
      \vskip1ex
      The \emph{Open} and \emph{Close} prices are the first and last trade prices recorded in the time bar.
      \vskip1ex
      The \emph{High} and \emph{Low} prices are the highest and lowest trade prices recorded in the time bar.
      \vskip1ex
      The \emph{Volume} is the total trading volume recorded in the time bar.
      \vskip1ex
      The \emph{OHLC} data format provides a way of efficiently compressing \emph{TAQ} data, while preserving information about price levels, volatility (range), and trading volumes.
      \vskip1ex
      In addition, evenly spaced \emph{OHLC} data allows for easier analysis of multiple time series, since the prices for different assets are given at the same moments in time.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,size="tiny",tidy=FALSE,eval=TRUE>>=
# Load package HighFreq
library(HighFreq)
head(HighFreq::SPY)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Plotting High Frequency \protect\emph{OHLC} Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Aggregating high frequency \emph{TAQ} data into \emph{OHLC} format with lower periodicity allows for data compression while maintaining some information about volatility.
      <<earl_ohlc_chart,echo=TRUE,eval=FALSE,fig.width=7,fig.height=6,fig.show='hide'>>=
# Load package HighFreq
library(HighFreq)
# Define symbol
symboln <- "SPY"
# Load OHLC data
dirout <- "/Users/jerzy/Develop/data/hfreq/scrub/"
symboln <- load(file.path(dirout, paste0(symboln, ".RData")))
interval <-"2013-11-11 09:30:00/2013-11-11 10:30:00"
chart_Series(SPY[interval], name=symboln)
      @
      The package \emph{HighFreq} contains both \emph{TAQ} data and \emph{Open-High-Low-Close} (\emph{OHLC}) data.
      \vskip1ex
      If you are not able to install package \emph{HighFreq} then download the file \texttt{hf\_data.RData} from the NYU share drive and load it.
    \column{0.5\textwidth}
      \includegraphics[width=0.45\paperwidth]{figure/earl_ohlc_chart-1}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Package \protect\emph{HighFreq} for Managing High Frequency Data}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package \emph{HighFreq} contains functions for managing high frequency time series data, such as:
      \begin{itemize}
        \item converting \emph{TAQ} data to \emph{OHLC} format,
        \item chaining and joining time series,
        \item scrubbing bad data,
        \item managing time zones and alligning time indices,
        \item aggregating data to lower frequency (periodicity),
        \item calculating rolling aggregations (VWAP, Hurst exponent, etc.),
        \item calculating seasonality aggregations,
        \item estimating volatility, skewness, and higher moments,
      \end{itemize}
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Install package HighFreq from github
devtools::install_github(repo="algoquant/HighFreq")
# Load package HighFreq
library(HighFreq)
# Get documentation for package HighFreq
# Get short description
packageDescription(HighFreq)
# Load help page
help(package=HighFreq)
# List all datasets in HighFreq
data(package=HighFreq)
# List all objects in HighFreq
ls("package:HighFreq")
# Remove HighFreq from search path
detach("package:HighFreq")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Datasets in Package \protect\emph{HighFreq}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package \emph{HighFreq} contains several high frequency time series, in \emph{xts} format, stored in a file called \texttt{hf\_data.RData}:
      \begin{itemize}
        \item a time series called \texttt{SPY\_TAQ}, containing a single day of \emph{TAQ} data for the \emph{SPY} ETF.
        \item three time series called \texttt{SPY}, \texttt{TLT}, and \texttt{VXX}, containing intraday \texttt{1}-minute \emph{OHLC} price bars for the \emph{SPY}, \emph{TLT}, and \emph{VXX} ETFs.
      \end{itemize}
      Even after the \emph{HighFreq} package is loaded, its datasets aren't loaded into the workspace, so they aren't listed in the workspace.
      \vskip1ex
      That's because the datasets in package \emph{HighFreq} are set up for \emph{lazy loading}, which means they can be called as if they were loaded, even though they're not loaded into the workspace.
      \vskip1ex
      The datasets in package \emph{HighFreq} can be loaded into the workspace using the function \texttt{data()}.
      \vskip1ex
      The data is set up for \emph{lazy loading}, so it doesn't require calling \texttt{data(hf\_data)} to load it into the workspace before calling it.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Load package HighFreq
library(HighFreq)
# You can see SPY when listing objects in HighFreq
ls("package:HighFreq")
# You can see SPY when listing datasets in HighFreq
data(package=HighFreq)
# But the SPY dataset isn't listed in the workspace
ls()
# HighFreq datasets are lazy loaded and available when needed
head(HighFreq::SPY)
# Load all the datasets in package HighFreq
data(hf_data)
# HighFreq datasets are now loaded and in the workspace
head(HighFreq::SPY)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Distribution of High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      High frequency returns exhibit \emph{large negative skewness} and \emph{very large kurtosis} (leptokurtosis), or fat tails.
      \vskip1ex
      Student's \emph{t-distribution} has fat tails, so it fits high frequency returns much better than the normal distribution.
      \vskip1ex
      The function \texttt{fitdistr()} from package \emph{MASS} fits a univariate distribution into a sample of data, by performing \emph{maximum likelihood} optimization.
      \vskip1ex
      The function \texttt{hist()} calculates and plots a histogram, and returns its data \emph{invisibly}.
      <<echo=(-(1:1)),eval=FALSE>>=
library(rutils)  # Load package rutils
# Calculate SPY percentage returns
ohlc <- HighFreq::SPY
nrows <- NROW(ohlc)
closep <- log(quantmod::Cl(ohlc))
retp <- rutils::diffit(closep)
colnames(retp) <- "SPY"
# Standardize raw returns to make later comparisons
retp <- (retp - mean(retp))/sd(retp)
# Calculate moments and perform normality test
sapply(c(var=2, skew=3, kurt=4), function(x) sum(retp^x)/nrows)
tseries::jarque.bera.test(retp)
# Fit SPY returns using MASS::fitdistr()
optiml <- MASS::fitdistr(retp, densfun="t", df=2)
loc <- optiml$estimate[1]
scalev <- optiml$estimate[2]
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_dist.png}
      <<echo=(-(1:2)),eval=FALSE>>=
x11(width=6, height=5)
par(mar=c(3, 3, 2, 1), oma=c(1, 1, 1, 1))
# Plot histogram of SPY returns
histp <- hist(retp, col="lightgrey", mgp=c(2, 1, 0),
  xlab="returns (standardized)", ylab="frequency", xlim=c(-3, 3),
  breaks=1e3, freq=FALSE, main="Distribution of High Frequency SPY Returns")
# lines(density(retp, bw=0.2), lwd=3, col="blue")
# Plot t-distribution function
curve(expr=dt((x-loc)/scalev, df=2)/scalev,
      type="l", lwd=3, col="red", add=TRUE)
# Plot the Normal probability distribution
curve(expr=dnorm(x, mean=mean(retp),
  sd=sd(retp)), add=TRUE, lwd=3, col="blue")
# Add legend
legend("topright", inset=0.05, bty="n",
  leg=c("t-distr", "normal"), y.intersp=0.1,
  lwd=6, lty=1, col=c("red", "blue"))
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Distribution of Aggregated High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The distribution of returns depends on the sampling frequency.
      \vskip1ex
      High frequency returns aggregated to a lower periodicity become less negatively skewed and less fat tailed, and closer to the normal distribution.
      \vskip1ex
      The function \texttt{xts::to.period()} converts a time series to a lower periodicity (for example from hourly to daily periodicity).
      <<echo=TRUE,eval=FALSE>>=
# Hourly SPY percentage returns
closep <- log(Cl(xts::to.period(x=ohlc, period="hours")))
retsh <- rutils::diffit(closep)
retsh <- (retsh - mean(retsh))/sd(retsh)
# Daily SPY percentage returns
closep <- log(Cl(xts::to.period(x=ohlc, period="days")))
retd <- rutils::diffit(closep)
retd <- (retd - mean(retd))/sd(retd)
# Calculate moments
sapply(list(minutely=retp, hourly=retsh, daily=retd),
       function(rets) {sapply(c(var=2, skew=3, kurt=4),
                function(x) mean(rets^x))
})  # end sapply
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_dist_agg.png}
      <<echo=(-(1:2)),eval=FALSE>>=
x11(width=6, height=5)
par(mar=c(3, 3, 2, 1), oma=c(1, 1, 1, 1))
# Plot densities of SPY returns
plot(density(retp, bw=0.4), xlim=c(-3, 3),
     lwd=3, mgp=c(2, 1, 0), col="blue",
     xlab="returns (standardized)", ylab="frequency",
     main="Density of High Frequency SPY Returns")
lines(density(retsh, bw=0.4), lwd=3, col="green")
lines(density(retd, bw=0.4), lwd=3, col="red")
# Add legend
legend("topright", inset=0.05, bty="n",
  leg=c("minutely", "hourly", "daily"), y.intersp=0.1,
  lwd=6, lty=1, col=c("blue", "green", "red"))
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Estimating Rolling Volatility of High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The volatility of high frequency returns can be inflated by large overnight returns.
      \vskip1ex
      The large overnight returns can be scaled down by dividing them by the overnight time interval.
        <<echo=TRUE,eval=FALSE>>=
# Calculate rolling volatility of SPY returns
ret2013 <- retp["2013-11-11/2013-11-15"]
# Calculate rolling volatility
lookb <- 11 # Look-back interval
endd <- seq_along(ret2013)
startp <- c(rep_len(1, lookb),
  endd[1:(NROW(endd)-lookb)])
endd[endd < lookb] <- lookb
vol_rolling <- sapply(seq_along(endd),
  function(it) sd(ret2013[startp[it]:endd[it]]))
vol_rolling <- xts::xts(vol_rolling, zoo::index(ret2013))
# Extract time intervals of SPY returns
indeks <- c(60, diff(xts::.index(ret2013)))
head(indeks)
table(indeks)
# Scale SPY returns by time intervals
ret2013 <- 60*ret2013/indeks
# Calculate scaled rolling volatility
vol_scaled <- sapply(seq_along(endd),
  function(it) sd(ret2013[startp[it]:endd[it]]))
vol_rolling <- cbind(vol_rolling, vol_scaled)
vol_rolling <- na.omit(vol_rolling)
sum(is.na(vol_rolling))
sapply(vol_rolling, range)
      @
    \column{0.5\textwidth}
      % \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_vol_rolling.png}
        <<echo=TRUE,eval=FALSE>>=
# Plot rolling volatility
x11(width=6, height=5)
plot_theme <- chart_theme()
plot_theme$col$line.col <- c("blue", "red")
chart_Series(vol_rolling, theme=plot_theme,
             name="Rolling Volatility with Overnight Spikes")
legend("topright", legend=colnames(vol_rolling),
  inset=0.1, bg="white", lty=1, lwd=6, y.intersp=0.1,
  col=plot_theme$col$line.col, bty="n")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Daily Volume and Volatility of High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Trading volumes typically rise together with market price volatility.
      \vskip1ex
      The function \texttt{apply.daily()} from package \texttt{xts} applies functions to time series over daily periods.
      \vskip1ex
      The function \texttt{calc\_var\_ohlc()} from package \texttt{HighFreq} calculates the variance of an \emph{OHLC} time series using range estimators.
      <<echo=TRUE,eval=FALSE>>=
# Volatility of SPY
sqrt(HighFreq::calcvar_ohlc(ohlc))
# Daily SPY volatility and volume
volatd <- sqrt(xts::apply.daily(ohlc, FUN=calcvar_ohlc))
colnames(volatd) <- ("SPY_volatility")
volumv <- quantmod::Vo(ohlc)
volumd <- xts::apply.daily(volumv, FUN=sum)
colnames(volumd) <- ("SPY_volume")
# Plot SPY volatility and volume
datav <- cbind(volatd, volumd)["2008/2009"]
colv <- colnames(datav)
dygraphs::dygraph(datav,
  main="SPY Daily Volatility and Trading Volume") %>%
  dyAxis("y", label=colv[1], independentTicks=TRUE) %>%
  dyAxis("y2", label=colv[2], independentTicks=TRUE) %>%
  dySeries(name=colv[1], axis="y", col="red", strokeWidth=3) %>%
  dySeries(name=colv[2], axis="y2", col="blue", strokeWidth=3)
      @
    \column{0.5\textwidth}
      % \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_daily_volat_volume.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Beta of Volume vs Volatility of High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      As a general empirical rule, the \emph{trading volume} $\upsilon$ in a given time period is roughly proportional to the \emph{volatility} of the returns $\sigma$: $\upsilon \propto \sigma$.
      \vskip1ex
      The regression of the \emph{log trading volume} versus the \emph{log volatility} fails the \emph{Durbin-Watson test} for the autocorrelation of residuals.
      \vskip1ex
      But the regression of the \emph{differences} passes the \emph{Durbin-Watson test}.
      <<echo=TRUE,eval=FALSE>>=
# Regress log of daily volume vs volatility
datav <- log(cbind(volumd, volatd))
colv <- colnames(datav)
dframe <- as.data.frame(datav)
formulav <- as.formula(paste(colv, collapse="~"))
regmod <- lm(formulav, data=dframe)
# Durbin-Watson test for autocorrelation of residuals
lmtest::dwtest(regmod)
# Regress diff log of daily volume vs volatility
dframe <- as.data.frame(rutils::diffit(datav))
regmod <- lm(formulav, data=dframe)
lmtest::dwtest(regmod)
summary(regmod)
plot(formulav, data=dframe, main="SPY Daily Trading Volume vs Volatility (log scale)")
abline(regmod, lwd=3, col="red")
mtext(paste("beta =", round(coef(regmod)[2], 3)), cex=1.2, lwd=3, side=2, las=2, adj=(-0.5), padj=(-7))
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_daily_volume_volat_reg.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Beta of Hourly Trading Volume vs Volatility}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Hourly aggregations of high frequency data also support the rule that the \emph{trading volume} is roughly proportional to the \emph{volatility} of the returns: $\upsilon \propto \sigma$.
      <<echo=TRUE,eval=FALSE>>=
# 60 minutes of data in lookb interval
lookb <- 60 # Look-back interval
vol2013 <- volumv["2013"]
ret2013 <- retp["2013"]
# Define end points with beginning stub
nrows <- NROW(ret2013)
nagg <- nrows %/% lookb
endd <- nrows-lookb*nagg + (0:nagg)*lookb
startp <- c(1, endd[1:(NROW(endd)-1)])
# Calculate SPY volatility and volume
datav <- sapply(seq_along(endd), function(it) {
  endp <- startp[it]:endd[it]
  c(volume=sum(vol2013[endp]),
    volatility=sd(ret2013[endp]))
})  # end sapply
datav <- t(datav)
datav <- rutils::diffit(log(datav))
dframe <- as.data.frame(datav)
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_hourly_volume_volat_reg.png}
      <<echo=TRUE,eval=FALSE>>=
formulav <- as.formula(paste(colnames(datav), collapse="~"))
regmod <- lm(formulav, data=dframe)
lmtest::dwtest(regmod)
summary(regmod)
plot(formulav, data=dframe,
     main="SPY Hourly Trading Volume vs Volatility (log scale)")
abline(regmod, lwd=3, col="red")
mtext(paste("beta =", round(coef(regmod)[2], 3)), cex=1.2, lwd=3, side=2, las=2, adj=(-0.5), padj=(-7))
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{High Frequency Returns in Trading Time}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{trading time} (volume clock) is the time measured by the level of \emph{trading volume}, with the \emph{volume clock} running faster in periods of higher \emph{trading volume}.
      \vskip1ex
      The time-dependent volatility of high frequency returns (\emph{heteroskedasticity}) produces their \emph{leptokurtosis} (large kurtosis, or fat tails).
      \vskip1ex
      The returns can be divided by the \emph{square root of the trading volumes} to obtain scaled returns over equal trading volumes.
      \vskip1ex
      But the returns should not be divided by very small volumes below a certain threshold.
      \vskip1ex
      The scaled returns have a smaller \emph{skewness} and \emph{kurtosis}, and they also have even higher autocorrelations than unscaled returns.
      <<echo=TRUE,eval=FALSE>>=
# Scale returns using volume (volume clock)
retsc <- ifelse(volumv > 1e4, retp/sqrt(volumv), 0)
retsc <- retsc/sd(retsc)
# Calculate moments of scaled returns
nrows <- NROW(retp)
sapply(list(retp=retp, retsc=retsc),
  function(rets) {sapply(c(skew=3, kurt=4),
           function(x) sum((rets/sd(rets))^x)/nrows)
})  # end sapply
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_scaled.png}
      <<echo=TRUE,eval=FALSE>>=
x11(width=6, height=5)
par(mar=c(3, 3, 2, 1), oma=c(1, 1, 1, 1))
# Plot densities of SPY returns
plot(density(retp), xlim=c(-3, 3),
     lwd=3, mgp=c(2, 1, 0), col="blue",
     xlab="returns (standardized)", ylab="frequency",
     main="Density of Volume-scaled High Frequency SPY Returns")
lines(density(retsc, bw=0.4), lwd=3, col="red")
curve(expr=dnorm, add=TRUE, lwd=3, col="green")
# Add legend
legend("topright", inset=0.05, bty="n", y.intersp=0.1,
  leg=c("minutely", "scaled", "normal"),
  lwd=6, lty=1, col=c("blue", "red", "green"))
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Autocorrelations of High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{Ljung-Box} test, tests if the autocorrelations of a time series are \emph{statistically significant}.
      \vskip1ex
      The \emph{null hypothesis} of the \emph{Ljung-Box} test is that the autocorrelations are equal to zero.
      \vskip1ex
      The \emph{Ljung-Box} statistic is small for time series that have \emph{statistically insignificant} autocorrelations.
      \vskip1ex
      The function \texttt{Box.test()} calculates the \emph{Ljung-Box} test and returns the test statistic and its p-value.
      \vskip1ex
      For \emph{minutely SPY} returns, the \emph{Ljung-Box} statistic is large and its \emph{p}-value is very small, so we can conclude that \emph{minutely SPY} returns have statistically significant autocorrelations.
      \vskip1ex
      For \emph{scaled minutely SPY} returns, the \emph{Ljung-Box} statistic is even larger, so its autocorrelations are even more statistically significant.
      \vskip1ex
      \emph{SPY} returns aggregated to longer time intervals are less autocorrelated.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Ljung-Box test for minutely SPY returns
Box.test(retp, lag=10, type="Ljung")
# Ljung-Box test for daily SPY returns
Box.test(retd, lag=10, type="Ljung")
# Ljung-Box test statistics for scaled SPY returns
sapply(list(retp=retp, retsc=retsc),
  function(rets) {
    Box.test(rets, lag=10, type="Ljung")$statistic
})  # end sapply
# Ljung-Box test statistics for aggregated SPY returns
sapply(list(minutely=retp, hourly=retsh, daily=retd),
  function(rets) {
    Box.test(rets, lag=10, type="Ljung")$statistic
})  # end sapply
      @
      The level of the autocorrelations depends on the sampling frequency, with higher frequency returns having more significant negative autocorrelations.
      \vskip1ex
      As the returns are aggregated to a lower periodicity, they become less autocorrelated, with daily returns having almost insignificant autocorrelations.
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Partial Autocorrelations of High Frequency Returns}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      High frequency minutely \emph{SPY} returns have statistically significant negative autocorrelations.
      \vskip1ex
      \emph{SPY} returns \emph{scaled} by the trading volumes have even more significant negative autocorrelations.
      <<echo=TRUE,eval=FALSE>>=
# Set plot parameters
x11(width=6, height=8)
par(mar=c(4, 4, 2, 1), oma=c(0, 0, 0, 0))
layout(matrix(c(1, 2), ncol=1), widths=c(6, 6), heights=c(4, 4))
# Plot the partial autocorrelations of minutely SPY returns
pacfl <- pacf(as.numeric(retp), lag=10,
     xlab="lag", ylab="partial autocorrelation", main="")
title("Partial Autocorrelations of Minutely SPY Returns", line=1)
# Plot the partial autocorrelations of scaled SPY returns
pacfs <- pacf(as.numeric(retsc), lag=10,
     xlab="lag", ylab="partial autocorrelation", main="")
title("Partial Autocorrelations of Scaled SPY Returns", line=1)
# Calculate the sums of partial autocorrelations
sum(pacfl$acf)
sum(pacfs$acf)
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_pacf.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Market Liquidity, Trading Volume and Volatility}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Market illiquidity is defined as the market price impact resulting from supply-demand imbalance.
      \vskip1ex
      Market liquidity $\mathcal{L}$ is proportional to the square root of the \emph{trading volume} $\upsilon$ divided by the price volatility $\sigma$:
      \begin{displaymath}
        \mathcal{L} \sim \frac{\sqrt{\upsilon}}{\sigma}
      \end{displaymath}
      Market illiquidity spiked during the May 6, 2010 \emph{flash crash}.
      \vskip1ex
      Research suggests that market crashes are caused by declining market liquidity:\\
      \href{http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2583743}{\emph{Donier et al., Why Do Markets Crash?}}
      <<echo=TRUE,eval=FALSE>>=
# Calculate market illiquidity
liquidv <- sqrt(volumd)/volatd
# Plot market illiquidity
x11(width=6, height=7) ; par(mfrow=c(2, 1))
plot_theme <- chart_theme()
plot_theme$col$line.col <- c("blue")
chart_Series(liquidv["2010"], theme=plot_theme,
  name="SPY Liquidity in 2010", plot=FALSE)
plot_theme$col$line.col <- c("red")
chart_Series(volatd["2010"],
  theme=plot_theme, name="SPY Volatility in 2010")
@
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_liquidity.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Intraday Seasonality of Volume and Volatility}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The volatility and trading volumes are typically higher at the beginning and end of the trading sessions.
      <<echo=TRUE,eval=FALSE>>=
# Calculate intraday time index with hours and minutes
datev <- format(zoo::index(retp), "%H:%M")
# Aggregate the mean volume
volumagg <- tapply(X=volumv, INDEX=datev, FUN=mean)
volumagg <- drop(volumagg)
# Aggregate the mean volatility
volagg <- tapply(X=retp^2, INDEX=datev, FUN=mean)
volagg <- sqrt(drop(volagg))
# Coerce to xts
datev <- as.POSIXct(paste(Sys.Date(), names(volumagg)))
volumagg <- xts::xts(volumagg, datev)
volagg <- xts::xts(volagg, datev)
# Plot seasonality of volume and volatility
x11(width=6, height=7) ; par(mfrow=c(2, 1))
plot_theme <- chart_theme()
plot_theme$col$line.col <- c("blue")
chart_Series(volumagg[c(-1, -NROW(volumagg))], theme=plot_theme,
  name="Intraday Seasonality of SPY Volume", plot=FALSE)
plot_theme$col$line.col <- c("red")
chart_Series(volagg[c(-1, -NROW(volagg))], theme=plot_theme,
  name="Intraday Seasonality of SPY Volatility")
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_intraday_volume_volat.png}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Intraday Seasonality of Liquidity and Volatility}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Market liquidity is typically the highest at the end of the trading session, and the lowest at the beginning.
      \vskip1ex
      The end of day spike in trading volumes and liquidity is driven by computer-driven investors liquidating their positions.
      <<echo=TRUE,eval=FALSE>>=
# Calculate market liquidity
liquidv <- sqrt(volumagg)/volagg
# Plot intraday seasonality of market liquidity
x11(width=6, height=7) ; par(mfrow=c(2, 1))
plot_theme <- chart_theme()
plot_theme$col$line.col <- c("blue")
chart_Series(liquidv[c(-1, -NROW(liquidv))], theme=plot_theme,
  name="Intraday Seasonality of SPY Liquidity", plot=FALSE)
plot_theme$col$line.col <- c("red")
chart_Series(volagg[c(-1, -NROW(volagg))], theme=plot_theme,
  name="Intraday Seasonality of SPY Volatility")
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{figure/hf_intraday_liquid_volat.png}
  \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Package \protect\emph{Rcpp} for Running \texttt{C++} Programs}


%%%%%%%%%%%%%%%
\subsection{Package \protect\emph{Rcpp} for Calling \texttt{C++} Programs from \texttt{R}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package \emph{Rcpp} allows calling \texttt{C++} functions from \texttt{R}, by compiling the \texttt{C++} code and creating \texttt{R} functions.
      \vskip1ex
      \emph{Rcpp} functions are \texttt{R} functions that were compiled from \texttt{C++} code using package \emph{Rcpp}.
      \vskip1ex
      \emph{Rcpp} functions are much faster than code written in \texttt{R}, so they're suitable for large numerical calculations.
      \vskip1ex
      The package \emph{Rcpp} relies on \emph{Rtools} for compiling the \texttt{C++} code: \\
      \hskip1em\url{https://cran.r-project.org/bin/windows/Rtools/}
      \vskip1ex
      You can learn more about the package \emph{Rcpp} here: \\
      \hskip1em\url{http://adv-r.had.co.nz/Rcpp.html}\\
      \hskip1em\url{http://www.rcpp.org/}\\
      \hskip1em\url{http://gallery.rcpp.org/}
    \column{0.5\textwidth}
      \vspace{-1em}
      \includegraphics[width=0.45\paperwidth]{image/bart_simpson_cpp.jpg}
      <<echo=TRUE,eval=FALSE>>=
# Verify that Rtools or XCode are working properly:
devtools::find_rtools()  # Under Windows
devtools::has_devel()
# Install the packages Rcpp and RcppArmadillo
install.packages(c("Rcpp", "RcppArmadillo"))
# Load package Rcpp
library(Rcpp)
# Get documentation for package Rcpp
# Get short description
packageDescription("Rcpp")
# Load help page
help(package="Rcpp")
# List all datasets in "Rcpp"
data(package="Rcpp")
# List all objects in "Rcpp"
ls("package:Rcpp")
# Remove Rcpp from search path
detach("package:Rcpp")
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Function \texttt{cppFunction()} for Compiling \texttt{C++} code}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The function \texttt{cppFunction()} compiles \texttt{C++} code into an \texttt{R} function.
      \vskip1ex
      The function \texttt{cppFunction()} creates an \texttt{R} function only for the current \texttt{R} session, and it must be recompiled for every new \texttt{R} session.
      \vskip1ex
      The function \texttt{sourceCpp()} compiles \texttt{C++} code contained in a file into \texttt{R} functions.
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Define Rcpp function
Rcpp::cppFunction("
  int times_two(int x)
    { return 2 * x;}
  ")  # end cppFunction
# Run Rcpp function
times_two(3)
# Source Rcpp functions from file
Rcpp::sourceCpp(file="/Users/jerzy/Develop/lecture_slides/scripts/mult_rcpp.cpp")
# Multiply two numbers
mult_rcpp(2, 3)
mult_rcpp(1:3, 6:4)
# Multiply two vectors
mult_vec_rcpp(2, 3)
mult_vec_rcpp(1:3, 6:4)
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Performing Loops in \protect\emph{Rcpp Sugar}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Loops written in \emph{Rcpp} can be two orders of magnitude faster than loops in \texttt{R}!
      \vskip1ex
      \emph{Rcpp Sugar} allows using \texttt{R}-style vectorized syntax in \emph{Rcpp} code.
      <<echo=TRUE,eval=FALSE>>=
# Define Rcpp function with loop
Rcpp::cppFunction("
double inner_mult(NumericVector x, NumericVector y) {
int xsize = x.size();
int ysize = y.size();
if (xsize != ysize) {
    return 0;
  } else {
    double total = 0;
    for(int i = 0; i < xsize; ++i) {
      total += x[i] * y[i];
  }
  return total;
  }
}")  # end cppFunction
# Run Rcpp function
inner_mult(1:3, 6:4)
inner_mult(1:3, 6:3)
# Define Rcpp Sugar function with loop
Rcpp::cppFunction("
double inner_sugar(NumericVector x, NumericVector y) {
  return sum(x * y);
}")  # end cppFunction
# Run Rcpp Sugar function
inner_sugar(1:3, 6:4)
inner_sugar(1:3, 6:3)
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Define R function with loop
inner_multr <- function(x, y) {
    sumv <- 0
    for(i in 1:NROW(x)) {
      sumv <- sumv + x[i] * y[i]
    }
    sumv
}  # end inner_multr
# Run R function
inner_multr(1:3, 6:4)
inner_multr(1:3, 6:3)
# Compare speed of Rcpp and R
library(microbenchmark)
summary(microbenchmark(
  rcode=inner_multr(1:10000, 1:10000),
  innerp=1:10000 %*% 1:10000,
  Rcpp=inner_mult(1:10000, 1:10000),
  sugar=inner_sugar(1:10000, 1:10000),
  times=10))[, c(1, 4, 5)]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Simulating Ornstein-Uhlenbeck Process Using \protect\emph{Rcpp}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      Simulating the Ornstein-Uhlenbeck Process in \emph{Rcpp} is about 30 times faster than in \texttt{R}!
      <<echo=TRUE,eval=FALSE>>=
# Define Ornstein-Uhlenbeck function in R
sim_our <- function(nrows=1000, priceq=5.0,
                    volat=0.01, theta=0.01) {
  retp <- numeric(nrows)
  pricev <- numeric(nrows)
  pricev[1] <- priceq
  for (i in 2:nrows) {
    retp[i] <- theta*(priceq - pricev[i-1]) + volat*rnorm(1)
    pricev[i] <- pricev[i-1] + retp[i]
  }  # end for
  pricev
}  # end sim_our
# Simulate Ornstein-Uhlenbeck process in R
priceq <- 5.0; sigmav <- 0.01
thetav <- 0.01; nrows <- 1000
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")  # Reset random numbers
ousim <- sim_our(nrows, priceq=priceq, volat=sigmav, theta=thetav)
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Define Ornstein-Uhlenbeck function in Rcpp
Rcpp::cppFunction("
NumericVector sim_oucpp(double priceq, 
                        double volat, 
                        double thetav, 
                        NumericVector innov) {
  int nrows = innov.size();
  NumericVector pricev(nrows);
  NumericVector retv(nrows);
  pricev[0] = priceq;
  for (int it = 1; it < nrows; it++) {
    retv[it] = thetav*(priceq - pricev[it-1]) + volat*innov[it-1];
    pricev[it] = pricev[it-1] + retv[it];
  }  // end for
  return pricev;
}")  # end cppFunction
# Simulate Ornstein-Uhlenbeck process in Rcpp
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")  # Reset random numbers
oucpp <- sim_oucpp(priceq=priceq,
  volat=sigmav, theta=thetav, innov=rnorm(nrows))
all.equal(ousim, oucpp)
# Compare speed of Rcpp and R
library(microbenchmark)
summary(microbenchmark(
  rcode=sim_our(nrows, priceq=priceq, volat=sigmav, theta=thetav),
  Rcpp=sim_oucpp(priceq=priceq, volat=sigmav, theta=thetav, innov=rnorm(nrows)),
  times=10))[, c(1, 4, 5)]
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{\protect\emph{Rcpp Attributes}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{Rcpp attributes} are instructions for the \texttt{C++} compiler, embedded in the \emph{Rcpp} code as \texttt{C++} comments, and preceded by the \texttt{"//"} symbol.
      \vskip1ex
      The \texttt{Rcpp::depends} attribute specifies additional \texttt{C++} library dependencies.
      \vskip1ex
      The \texttt{Rcpp::export} attribute specifies that a function should be exported to \texttt{R}, where it can be called as an \texttt{R} function.
      \vskip1ex
      Only functions which are preceded by the \texttt{Rcpp::export} attribute are exported to \texttt{R}.
      \vskip1ex
      The function \texttt{sourceCpp()} compiles \texttt{C++} code contained in a file into \texttt{R} functions.
      <<echo=TRUE,eval=FALSE>>=
# Source Rcpp function for Ornstein-Uhlenbeck process from file
Rcpp::sourceCpp(file="/Users/jerzy/Develop/lecture_slides/scripts/sim_ou.cpp")
# Simulate Ornstein-Uhlenbeck process in Rcpp
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")  # Reset random numbers
oucpp <- sim_oucpp(priceq=priceq,
  volat=sigmav,
  theta=thetav,
  innov=rnorm(nrows))
all.equal(ousim, oucpp)
# Compare speed of Rcpp and R
library(microbenchmark)
summary(microbenchmark(
  rcode=sim_our(nrows, priceq=priceq, volat=sigmav, theta=thetav),
  Rcpp=sim_oucpp(priceq=priceq, volat=sigmav, theta=thetav, innov=rnorm(nrows)),
  times=10))[, c(1, 4, 5)]
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Rcpp header with information for C++ compiler
#include <Rcpp.h> // include Rcpp C++ header files
using namespace Rcpp; // use Rcpp C++ namespace

// The function sim_oucpp() simulates an Ornstein-Uhlenbeck process
// export the function roll_maxmin() to R
// [[Rcpp::export]]
NumericVector sim_oucpp(double priceq,
                          double volat,
                          double thetav,
                          NumericVector innov) {
  int(nrows = innov.size();
  NumericVector pricev*nrows);
  NumericVector retp*nrows);
  pricev[0] = priceq;
  for (int it = 1; it < nrows; it++) {
    retp[it] = thetav*(priceq - pricev[it-1]) + volat*innov[it-1];
    pricev[it] = pricev[it-1] + retp[it];
  }  // end for
  return pricev;
}  // end sim_oucpp
    \end{lstlisting}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Generating Random Numbers Using Logistic Map in \protect\emph{Rcpp}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The \emph{logistic map} in \emph{Rcpp} is about seven times faster than the loop in \texttt{R}, and even slightly faster than the standard \texttt{runif()} function in \texttt{R}!
      <<echo=TRUE,eval=FALSE>>=
# Calculate uniformly distributed pseudo-random sequence
unifun <- function(seedv, nrows=10) {
  datav <- numeric(nrows)
  datav[1] <- seedv
  for (i in 2:nrows) {
    datav[i] <- 4*datav[i-1]*(1-datav[i-1])
  }  # end for
  acos(1-2*datav)/pi
}  # end unifun

# Source Rcpp functions from file
Rcpp::sourceCpp(file="/Users/jerzy/Develop/lecture_slides/scripts/unifun.cpp")
# Microbenchmark Rcpp code
library(microbenchmark)
summary(microbenchmark(
  rcode=runif(1e5),
  rloop=unifun(0.3, 1e5),
  Rcpp=unifuncpp(0.3, 1e5),
  times=10))[, c(1, 4, 5)]
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Rcpp header with information for C++ compiler
#include <Rcpp.h> // include Rcpp C++ header files
using namespace Rcpp; // use Rcpp C++ namespace

// This is a simple example of exporting a C++ function to R.
// You can source this function into an R session using the
// function Rcpp::sourceCpp()
// (or via the Source button on the editor toolbar).
// Learn more about Rcpp at:
//
//   http://www.rcpp.org/
//   http://adv-r.had.co.nz/Rcpp.html
//   http://gallery.rcpp.org/

// function unifun() produces a vector of
// uniformly distributed pseudo-random numbers
// [[Rcpp::export]]
NumericVector unifuncpp(double seedv, int(nrows) {
// define pi
static const double pi = 3.14159265;
// allocate output vector
  NumericVector datav(nrows);
// initialize output vector
  datav[0] = seedv;
// perform loop
  for (int i=1; i < nrows; ++i) {
    datav[i] = 4*datav[i-1]*(1-datav[i-1]);
  }  // end for
// rescale output vector and return it
  return acos(1-2*datav)/pi;
}
    \end{lstlisting}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Package \protect\emph{RcppArmadillo} for Fast Linear Algebra}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package \emph{RcppArmadillo} allows calling from \texttt{R} the high-level \emph{Armadillo} \texttt{C++} linear algebra library.
      \vskip1ex
      \emph{Armadillo} provides ease of use and speed, with syntax similar to \emph{Matlab}.
      \vskip1ex
      \emph{RcppArmadillo} functions are often faster than even compiled \texttt{R} functions, because they use better optimized \texttt{C++} code:\\
      \url{http://arma.sourceforge.net/speed.html}\\
      \vskip1ex
      You can learn more about \emph{RcppArmadillo}: \\
      \tiny \url{http://arma.sourceforge.net/}\\
      \tiny \url{http://dirk.eddelbuettel.com/code/rcpp.armadillo.html}\\
      \tiny \url{https://cran.r-project.org/web/packages/\emph{RcppArmadillo}/index.html}\\
      \tiny \url{https://github.com/RcppCore/\emph{RcppArmadillo}}
      <<echo=TRUE,eval=FALSE>>=
library(RcppArmadillo)
# Source Rcpp functions from file
Rcpp::sourceCpp(file="/Users/jerzy/Develop/lecture_slides/scripts/armadillo_functions.cpp")
vec1 <- runif(1e5)
vec2 <- runif(1e5)
inner_vec(vec1, vec2)
vec1 %*% vec2
      @
    \column{0.5\textwidth}
      \vspace{-2em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Rcpp header with information for C++ compiler
#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;
// [[Rcpp::depends(RcppArmadillo)]]

// The function inner_vec() calculates the inner (dot) product of two vectors.
// It uses \emph{RcppArmadillo}.
//' @export
// [[Rcpp::export]]
double inner_vec(arma::vec vec1, arma::vec vec2) {
  return arma::dot(vec1, vec2);
}  // end inner_vec

// The function inner_mat() calculates the inner (dot) product of a matrix
// with two vectors.
// It accepts pointers to the matrix and vectors, and returns a double.
// It uses \emph{RcppArmadillo}.
//' @export
// [[Rcpp::export]]
double inner_mat(const arma::vec& vecv2, const arma::mat& matv, const arma::vec& vecv1) {
  return arma::as_scalar(trans(vecv2) * (matv * vecv1));
}  // end inner_mat
    \end{lstlisting}
      \vspace{-1.5em}
      <<echo=TRUE,eval=FALSE>>=
# Microbenchmark \emph{RcppArmadillo} code
summary(microbenchmark(
  rcpp = inner_vec(vec1, vec2),
  rcode = (vec1 %*% vec2),
  times=100))[, c(1, 4, 5)]  # end microbenchmark summary
# Microbenchmark shows:
# inner_vec() is several times faster than %*%, especially for longer vectors.
#     expr     mean   median
# 1 inner_vec 110.7067 110.4530
# 2 rcode 585.5127 591.3575
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Simulating \protect\emph{ARIMA} Processes Using \protect\emph{RcppArmadillo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{ARIMA} processes can be simulated using \emph{RcppArmadillo} even faster than by using the function \texttt{filter()}.
      <<echo=TRUE,eval=FALSE>>=
# Source Rcpp functions from file
Rcpp::sourceCpp(file="/Users/jerzy/Develop/lecture_slides/scripts/sim_arima.cpp")
# Define AR(2) coefficients
coeff <- c(0.9, 0.09)
nrows <- 1e4
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
innov <- rnorm(nrows)
# Simulate ARIMA using filter()
arimar <- filter(x=innov, filter=coeff, method="recursive")
# Simulate ARIMA using sim_ar()
innov <- matrix(innov)
coeff <- matrix(coeff)
arimav <- sim_ar(coeff, innov)
all.equal(drop(arimav), as.numeric(arimar))
# Microbenchmark \emph{RcppArmadillo} code
summary(microbenchmark(
  rcpp = sim_ar(coeff, innov),
  filter = filter(x=innov, filter=coeff, method="recursive"),
  times=100))[, c(1, 4, 5)]  # end microbenchmark summary
      @
    \column{0.5\textwidth}
      % \vspace{-2em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Rcpp header with information for C++ compiler
#include <RcppArmadillo.h> // include C++ header file from Armadillo library
using namespace arma; // use C++ namespace from Armadillo library
// declare dependency on RcppArmadillo
// [[Rcpp::depends(RcppArmadillo)]]

//' @export
// [[Rcpp::export]]
arma::vec sim_ar(const arma::vec& innov, const arma::vec& coeff) {
  uword nrows = innov.n_elem;
  uword lookb = coeff.n_elem;
  arma::vec arimav[nrows);

  // startup period
  arimav(0) = innov(0);
  arimav(1) = innov(1) + coeff(lookb-1) * arimav(0);
  for (uword it = 2; it < lookb-1; it++) {
    arimav(it) = innov(it) + arma::dot(coeff.subvec(lookb-it, lookb-1), arimav.subvec(0, it-1));
  }  // end for

  // remaining periods
  for (uword it = lookb; it < nrows; it++) {
    arimav(it) = innov(it) + arma::dot(coeff, arimav.subvec(it-lookb, it-1));
  }  // end for

  return arimav;
}  // end sim_arima
    \end{lstlisting}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Fast Matrix Algebra Using \protect\emph{RcppArmadillo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{RcppArmadillo} functions can be made even faster by operating on pointers to matrices and performing calculations in place, without copying large matrices.
      \vskip1ex
      \emph{RcppArmadillo} functions can be compiled using the same \emph{Rtools} as those for \emph{Rcpp} functions:\\
      \hskip1em\url{https://cran.r-project.org/bin/windows/Rtools/}
      <<echo=TRUE,eval=FALSE>>=
library(RcppArmadillo)
# Source Rcpp functions from file
Rcpp::sourceCpp(file="/Users/jerzy/Develop/lecture_slides/scripts/armadillo_functions.cpp")
matv <- matrix(runif(1e5), nc=1e3)
# Center matrix columns using apply()
matd <- apply(matv, 2, function(x) (x-mean(x)))
# Center matrix columns in place using Rcpp demeanr()
demeanr(matv)
all.equal(matd, matv)
# Microbenchmark \emph{RcppArmadillo} code
library(microbenchmark)
summary(microbenchmark(
  rcode = (apply(matv, 2, mean)),
  rcpp = demeanr(matv),
  times=100))[, c(1, 4, 5)]  # end microbenchmark summary
# Perform matrix inversion
# Create random positive semi-definite matrix
matv <- matrix(runif(25), nc=5)
matv <- t(matv) %*% matv
# Invert the matrix
matrixinv <- solve(matv)
inv_mat(matv)
all.equal(matrixinv, matv)
# Microbenchmark \emph{RcppArmadillo} code
summary(microbenchmark(
  rcode = solve(matv),
  rcpp = inv_mat(matv),
  times=100))[, c(1, 4, 5)]  # end microbenchmark summary
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Rcpp header with information for C++ compiler
#include <RcppArmadillo.h> // include C++ header file from Armadillo library
using namespace arma; // use C++ namespace from Armadillo library
// declare dependency on RcppArmadillo
// [[Rcpp::depends(RcppArmadillo)]]

// Examples of \emph{RcppArmadillo} functions below

// The function demeanr() calculates a matrix with centered columns.
// It accepts a pointer to a matrix and operates on the matrix in place.
// It returns the number of columns of the input matrix.
// It uses \emph{RcppArmadillo}.
//' @export
// [[Rcpp::export]]
int demeanr(arma::mat& matv) {
  for (uword i = 0; i < matv.n_cols; i++) {
    matv.col(i) -= arma::mean(matv.col(i));
  }  // end for
  return matv.n_cols;
}  // end demeanr

// The function inv_mat() calculates the inverse of symmetric positive
// definite matrix.
// It accepts a pointer to a matrix and operates on the matrix in place.
// It returns the number of columns of the input matrix.
// It uses \emph{RcppArmadillo}.
//' @export
// [[Rcpp::export]]
double inv_mat(arma::mat& matv) {
  matv = arma::inv_sympd(matv);
  return matv.n_cols;
}  // end inv_mat
    \end{lstlisting}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Fast Correlation Matrix Inverse Using \protect\emph{RcppArmadillo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      \emph{RcppArmadillo} can be used to quickly calculate the reduced inverse of correlation matrices.
      <<echo=TRUE,eval=FALSE>>=
library(RcppArmadillo)
# Source Rcpp functions from file
Rcpp::sourceCpp("/Users/jerzy/Develop/lecture_slides/scripts/HighFreq.cpp")
# Calculate matrix of random returns
matv <- matrix(rnorm(300), nc=5)
# Reduced inverse of correlation matrix
dimax <- 4
cormat <- cor(matv)
eigend <- eigen(cormat)
invmat <- eigend$vectors[, 1:dimax] %*%
  (t(eigend$vectors[, 1:dimax]) / eigend$values[1:dimax])
# Reduced inverse using \emph{RcppArmadillo}
invarma <- calc_inv(cormat, dimax=dimax)
all.equal(invmat, invarma)
# Microbenchmark \emph{RcppArmadillo} code
library(microbenchmark)
summary(microbenchmark(
  rcode = {eigend <- eigen(cormat)
      eigend$vectors[, 1:dimax] %*% (t(eigend$vectors[, 1:dimax]) / eigend$values[1:dimax])},
  rcpp = calc_inv(cormat, dimax=dimax),
  times=100))[, c(1, 4, 5)]  # end microbenchmark summary
      @
    \column{0.5\textwidth}
      \vspace{-1em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Rcpp header with information for C++ compiler
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>  // include Rcpp C++ header files
using namespace stdev;
using namespace Rcpp; // use Rcpp C++ namespace
using namespace arma;

//' @export
// [[Rcpp::export]]
arma::mat calc_inv(const arma::mat& matv, 
                   arma::uword dimax = 0, // Max number of PCA for dimension reduction
                   double eigen_thresh = 0.01) { // Threshold for discarding small singular values

  // Allocate SVD variables
  arma::vec svdval;  // Singular values
  arma::mat svdu, svdv;  // Singular matrices
  // Calculate the SVD
  arma::svd(svdu, svdval, svdv, tseries);
  // Calculate the number of non-small singular values
  arma::uword svdnum = arma::sum(svdval > eigen_thresh*arma::sum(svdval));
  
  // If no regularization then set dimax to (svdnum - 1)
  if (dimax == 0) {
    // Set dimax
    dimax = svdnum - 1;
  } else {
    // Adjust dimax
    dimax = stdev::min(dimax - 1, svdnum - 1);
  }  // end if
  
  // Remove all small singular values
  svdval = svdval.subvec(0, dimax);
  svdu = svdu.cols(0, dimax);
  svdv = svdv.cols(0, dimax);
  
  // Calculate the reduced inverse from the SVD decomposition
  return svdv*arma::diagmat(1/svdval)*svdu.t();
  
}  // end calc_inv
    \end{lstlisting}
    % \column{0.5\textwidth}
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Portfolio Optimization Using \protect\emph{RcppArmadillo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  % \begin{columns}[T]
    % \column{0.5\textwidth}
      Fast portfolio optimization using matrix algebra can be implemented using \protect\emph{RcppArmadillo}.
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
// Fast portfolio optimization using matrix algebra and \emph{RcppArmadillo}
arma::vec calc_weights(const arma::mat& returns, // Asset returns
                       Rcpp::List controlv) { // List of portfolio optimization parameters
  
  // Apply different calculation methods for weights
  switch(calc_method(method)) {
  case methodenum::maxsharpe: {
    // Mean returns of columns
    arma::vec colmeans = arma::trans(arma::mean(returns, 0));
    // Shrink colmeans to the mean of returns
    colmeans = ((1-alpha)*colmeans + alpha*arma::mean(colmeans));
    // Calculate weights using reduced inverse
    weights = calc_inv(covmat, dimax, eigen_thresh)*colmeans;
    break;
  }  // end maxsharpe
  case methodenum::maxsharpemed: {
    // Median returns of columns
    arma::vec colmeans = arma::trans(arma::median(returns, 0));
    // Shrink colmeans to the median of returns
    colmeans = ((1-alpha)*colmeans + alpha*arma::median(colmeans));
    // Calculate weights using reduced inverse
    weights = calc_inv(covmat, dimax, eigen_thresh)*colmeans;
    break;
  }  // end maxsharpemed
  case methodenum::minvarlin: {
    // Minimum variance weights under linear constraint
    // Multiply reduced inverse times unit vector
    weights = calc_inv(covmat, dimax, eigen_thresh)*arma::ones(ncols);
    break;
  }  // end minvarlin
  case methodenum::minvarquad: {
    // Minimum variance weights under quadratic constraint
    // Calculate highest order principal component
    arma::vec eigenval;
    arma::mat eigenvec;
    arma::eig_sym(eigenval, eigenvec, covmat);
    weights = eigenvec.col(ncols-1);
    break;
  }  // end minvarquad
  case methodenum::sharpem: {
    // Momentum weights equal to Sharpe ratios
    // Mean returns of columns
    arma::vec colmeans = arma::trans(arma::mean(returns, 0));
    // Standard deviation of columns
    arma::vec colsd = arma::sqrt(covmat.diag());
    colsd.replace(0, 1);
    // Momentum weights equal to Sharpe ratios
    weights = colmeans/colsd;
    break;
  }  // end sharpem
  case methodenum::kellym: {
    // Momentum weights equal to Kelly ratios
    // Mean returns of columns
    arma::vec colmeans = arma::trans(arma::mean(returns, 0));
    // Variance of columns
    arma::vec colvar = covmat.diag();
    colvar.replace(0, 1);
    // Momentum weights equal to Kelly ratios
    weights = colmeans/colvar;
    break;
  }  // end kellym
  case methodenum::robustm: {
    // Momentum weights equal to robust Sharpe ratios
    // Median returns of columns
    arma::vec colmeans = arma::trans(arma::median(returns, 0));
    // Standard deviation of columns
    arma::vec colsd = arma::sqrt(covmat.diag());
    colsd.replace(0, 1);
    // Momentum weights equal to robust Sharpe ratios
    colmeans = colmeans/colsd;
    break;
  }  // end robustm
  case methodenum::quantile: {
    // Momentum weights equal to sum of quantiles for columns
    arma::vec levels = {confl, 1-confl};
    weights = arma::conv_to<vec>::from(arma::sum(arma::quantile(returns, levels, 0), 0));
    break;
  }  // end quantile
  default : {
    cout << "Warning: Invalid method parameter: " << method << endl;
    return arma::ones(ncols);
  }  // end default
  }  // end switch
  
  if (rankw == TRUE) {
    // Convert the weights to their ranks
    weights = arma::conv_to<vec>::from(calc_ranks_stl(weights));
  }  // end if
  
  if (centerw == TRUE) {
    // Center the weights so their sum is equal to zero
    weights = (weightv - arma::mean(weights));
  }  // end if
  
  // Apply different scaling methods for weights
  switch(calc_method(scalew)) {
  case methodenum::voltarget: {
    // Scale the weights so the portfolio has the volatility equal to vol_target
    weights = weights*vol_target/arma::stddev(returns*weightv);
    break;
  }  // end voltarget
  case methodenum::voleqw: {
    // Scale the weights to the volatility of the equal weight portfolio
    weights = weights*arma::stddev(arma::mean(returns, 1))/arma::stddev(returns*weightv);
    break;
  }  // end voleqw
  case methodenum::sumone: {
    // Scale the weights so their sum of squares is equal to one
    weights = weights/arma::sum(weightv*arma::ones(ncols));
    break;
  }  // end sumone
  case methodenum::sumsq: {
    // Scale the weights so their sum of squares is equal to one
    weights = weights/stdev::sqrt(arma::sum(square(weights)));
    break;
  }  // end sumsq
  default : {
    // No scaling
    break;
  }  // end default
  }  // end switch
  
  return weights;
  
}  // end calc_weights
    \end{lstlisting}
    % \column{0.5\textwidth}
      % \vspace{-1em}
  % \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Strategy Backtesting Using \protect\emph{RcppArmadillo}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  % \begin{columns}[T]
    % \column{0.5\textwidth}
      Fast backtesting of strategies can be implemented using \protect\emph{RcppArmadillo}.
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
arma::mat back_test(const arma::mat& excess, // Asset excess returns
                    const arma::mat& returns, // Asset returns
                    Rcpp::List controlv, // List of portfolio optimization model parameters
                    arma::uvec startp, // Start points
                    arma::uvec endd, // End points
                    double lambdaf = 0.0, // Decay factor for averaging the portfolio weights
                    double coeff = 1.0, // Multiplier of strategy returns
                    double bidask = 0.0) { // The bid-ask spread
  
  double lambda1 = 1-lambdaf;
  arma::uword nweights = returns.n_cols;
  arma::vec weights(nweights, fill::zeros);
  arma::vec weights_past = ones(nweights)/stdev::sqrt(nweights);
  arma::mat pnls = zeros(returns.n_rows, 1);

  // Perform loop over the end points
  for (arma::uword it = 1; it < endd.size(); it++) {
    // cout << "it: " << it << endl;
    // Calculate the portfolio weights
    weights = coeff*calc_weights(excess.rows(startp(it-1), endd(it-1)), controlv);
    // Calculate the weights as the weighted sum with past weights
    weights = lambda1*weights + lambdaf*weights_past;
    // Calculate out-of-sample returns
    pnls.rows(endd(it-1)+1, endd(it)) = returns.rows(endd(it-1)+1, endd(it))*weights;
    // Add transaction costs
    pnls.row(endd(it-1)+1) -= bidask*sum(abs(weightv - weights_past))/2;
    // Copy the weights
    weights_past = weights;
  }  // end for
  
  // Return the strategy pnls
  return pnls;
  
}  // end back_test
    \end{lstlisting}
    % \column{0.5\textwidth}
      % \vspace{-1em}
  % \end{columns}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Package \protect\emph{reticulate} for Running \texttt{Python} Programs}


%%%%%%%%%%%%%%%
\subsection{Package \protect\emph{reticulate} for Running \texttt{Python} from \texttt{RStudio}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
  \begin{columns}[T]
    \column{0.5\textwidth}
      The package \emph{reticulate} allows running \texttt{Python} functions and scripts from \texttt{RStudio}.
      \vskip1ex
      The package \emph{reticulate} relies on \texttt{Python} for interpreting the \texttt{Python} code.
      \vskip1ex
      You must set your Global Options in \texttt{RStudio} to your \texttt{Python} executable, for example:\\
      /Library/Frameworks/Python.framework/Versions/3.10/bin/python3.10
      \vskip1ex
      You can learn more about the package \emph{reticulate} here: \\
      \hskip1em\url{https://rstudio.github.io/reticulate/}
    \column{0.5\textwidth}
      \vspace{-1em}
      <<echo=TRUE,eval=FALSE>>=
# Install package reticulate
install.packages("reticulate")
# Start Python session
reticulate::repl_python()
# Exit Python session
exit
      @
  \end{columns}
\end{block}

\end{frame}


%%%%%%%%%%%%%%%
\subsection{Running \texttt{Python} Under \protect\emph{reticulate}}
\begin{frame}[fragile,t]{\subsecname}
\vspace{-1em}
\begin{block}{}
      \vspace{-2em}
      \begin{lstlisting}[language=R,basicstyle=\tiny\ttfamily\bfseries,backgroundcolor=\color{anti_flashwhite},showstringspaces=FALSE]
"""
Script for loading OHLC data from a CSV file and plotting a candlestick plot.
"""
# Import packages 
import pandas as pd
import numpy as np
import plotly.graph_objects as go
# Load OHLC data from csv file - the time index is formatted inside read_csv()
symbol = "SPY"
range = "day"
filename = "/Users/jerzy/Develop/data/" + symbol + "_" + range + ".csv"
ohlc = pd.read_csv(filename)
datev = ohlc.Date
# Calculate log stock prices
ohlc[["Open", "High", "Low", "Close"]] = np.log(ohlc[["Open", "High", "Low", "Close"]])
# Calculate moving average
lookback = 55
closep = ohlc.Close
pricema = closep.ewm(span=lookback, adjust=False).mean()
# Plotly simple candlestick with moving average
# Create empty graph object
plotfig = go.Figure()
# Add trace for candlesticks
plotfig = plotfig.add_trace(go.Candlestick(x=datev,
  open=ohlc.Open, high=ohlc.High, low=ohlc.Low, close=ohlc.Close, 
  name=symbol+" Log OHLC Prices", showlegend=False))
# Add trace for moving average
plotfig = plotfig.add_trace(go.Scatter(x=datev, y=pricema, 
  name="Moving Average", line=dict(color="blue")))
# Customize plot
plotfig = plotfig.update_layout(title=symbol + " Log OHLC Prices", 
  title_font_size=24, title_font_color="blue", yaxis_title="Price", 
  font_color="black", font_size=18, xaxis_rangeslider_visible=False)
# Customize legend
plotfig = plotfig.update_layout(legend=dict(x=0.2, y=0.9, traceorder="normal", 
  itemsizing="constant", font=dict(family="sans-serif", size=18, color="blue")))
# Render the plot
plotfig.show()
    \end{lstlisting}
\end{block}

\end{frame}



%%%%%%%%%%%%%%%
\section{Homework Assignment}


%%%%%%%%%%%%%%%
\subsection{\secname}
\begin{frame}[t]{\subsecname}
\vspace{-1em}
\begin{block}{No homework!}
  \hskip10.0em\includegraphics[scale=0.1]{image/smile.png}
\end{block}

\end{frame}


\end{document}