MAthesis_Statistics_Berkeley_Spring2014_HenryWong.rnw

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% The Sweave file MAthesis_Statistics_Berkeley_Spring2014_HenryWong.rnw

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[masters]{ucbthesis}
\usepackage{biblatex}

% Double spacing, if you want it.
% \def\dsp{\def\baselinestretch{2.0}\large\normalsize}
% \dsp

% If the Grad. Division insists that the first paragraph of a section
% be indented (like the others), then include this line:
% \usepackage{indentfirst}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage[sc]{mathpazo}
\usepackage[T1]{fontenc}
\setcounter{secnumdepth}{5}
\setcounter{tocdepth}{2}
\usepackage{url}
\usepackage[unicode=true,pdfusetitle,
 bookmarks=true,bookmarksnumbered=true,bookmarksopen=true,bookmarksopenlevel=2,
 breaklinks=false,pdfborder={0 0 1},backref=false,colorlinks=false]
 {hyperref}
\hypersetup{
 pdfstartview={XYZ null null 1}}
\usepackage{breakurl}
\usepackage{float}
\floatstyle{boxed}
\restylefloat{figure}
\usepackage{amsmath}
\usepackage{url}
\usepackage{multirow}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newtheorem{theorem}{Jibberish}
\bibliography{references}
\hyphenation{mar-gin-al-ia}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

% Declarations for Front Matter

\title{Statistical Validity of the Fama French Carhart Four Model}
\author{Henry Ka Ming Wong}
\degreesemester{Spring}
\degreeyear{2014}
\degree{Master of Arts}
\cochairs{Professor Lisa R. Goldberg}{Professor Michael J. Klass}
\othermembers{Professor Robert M. Anderson}
\numberofmembers{3}
\prevdegrees{B.S. (University of California, Los Angeles) 2004 \\
  M.S. (University of California, Berkeley) 2013}
\field{Statistics}
\campus{Berkeley}

% For a masters thesis, uncomment (remove the % at the beginning of)
% the following line.  This affects the title and approval pages,
% which by default calls this a "dissertation", not a "thesis".

%\itsamasters

% The title page generated by LaTeX is now acceptable for handing in.
% (This was not always the case).

\maketitle
\approvalpage
\copyrightpage

\begin{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{dedication}
I dedicate this thesis to my immediate family members: my parents, my brothers David and Alvin, and my sister-in-law Lily. Thank you for all of their love and support throughout my life. In particular, I would like to thank my parents for always believing in me, their continuous love, and their support in all my life decisions. Without whom I could not have made it here. My bothers David and Alvin have never left my side and are very special.
\\

I also dedicate this thesis to my friends who have supported me throughout the process. I will always appreciate all that they have done, especially Ben Chan, Kit Chan, and their wonderful children Ada and Paul, for being supportive and encouraging as always. They have been my best cheerleaders.
\end{dedication}

\tableofcontents
\clearpage
\listoffigures
\clearpage
\listoftables

\begin{acknowledgements}
My first and sincere appreciation goes to Lisa R. Goldberg and Robert M. Anderson, my advisors, for all I have learned from them and for their continuous help and support in all stages of this thesis. I especially thank them for suggesting the research topic for this thesis, and for encouraging and helping me shape my interest and ideas. I appreciate all their contribution of time and ideas to make my research experience productive and stimulating. The joy and enthusiasm they have for their own research have been contagious and motivational for me.
\\

I would like to express my deep gratitude and respect to Michael J. Klass, whose advices, insight, and friendship have been invaluable to me. I very much appreciate his enthusiasm, intensity, and willingness to open to new ideas. I greatly value the close personal rapport that Michael and I have forged over the last two years. I quite simply cannot imagine a better advisor.
\\

In addition, I would like to thank Christopher Paciorek for his interesting "Statistical Computing" course, which advanced my skills in R, Unix, and $\LaTeX$ that I used for this research. I also want to thank his ongoing support in the software programs even after the course ended.
\\

In the end, I want to thank Victoria Stodden for her suggestion in using GitHub to facilitate reproducible computational results and share the research for this thesis. Her suggestion has helped me pave the path of sharing codes and research ideas with scholars in my field. 
\end{acknowledgements}

\end{frontmatter}

\pagestyle{headings}

% (Optional) \part{First Part}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<setup, include=FALSE, cache=FALSE>>=
library(knitr)
# set global chunk options
opts_chunk$set(fig.path='figure/fig-', fig.show='asis')
options(replace.assign=TRUE,width=90)
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run the R codes for the required functions and variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

<<echo=FALSE, message=FALSE, warning=FALSE>>=
#############################################################################
### The Beginning of R codes for functions and variables
#############################################################################
### Import the data ###
#############################################################################

library(sandwich)
library(lmtest)
library(tseries)
library(xtable)
library(nortest)
library(car)

### Import the data
portfolios.e <- read.table("10_Portfolios_AveEquWeight.txt", header=T,
quote="\"")
portfolios.v <- read.table("10_Portfolios_AveValWeight.txt", header=T,
quote="\"")
researchfactor <- read.table("F-F_Research_Data_Factors.txt", header=T,
quote="\"")
momfactor <- read.table("F-F_Momentum_Factor.TXT", header=T,
quote="\"")

### Create function to combine the data
dataFun <- function(portfolios, researchfactor, momfactor){
    
    ### Combine the data
    data <- cbind(portfolios, researchfactor[,2:4], momfactor[,2],
    researchfactor[,5])
    
    ### Rename the variables
    names(data)[2:11] <- paste("Portfolio", 1:10, sep="")
    names(data)[12] <- "RMRF"
    names(data)[15:16] <- c("MOM", "RF")
    
    ### Convert the data in percentage
    data[,2:16] <- data[,2:16] / 100
    
    ### Compute the excess portfolio returns
    data[,2:11] <- data[,2:11] - data[,16]
    
    ### Convert the month from numeric to date format
    month.temp <- as.character(data[,1])
    data[,1] <- paste(substring(month.temp, 1, 4),
    substring(month.temp, 5, 6), "01", sep="-")
    data[,1] <- as.Date(data[,1])
    
    return (data)
}

### Get the average equal weighted (AEW) data with excess return
data.e <- dataFun(portfolios.e, researchfactor, momfactor)

### Get the average value weighted (AVW) data with excess return
data.v <- dataFun(portfolios.v, researchfactor, momfactor)

### Data with real return
pdata.e <- data.e
pdata.v <- data.v
pdata.e[,2:11] <- pdata.e[,2:11] + pdata.e[,16]
pdata.v[,2:11] <- pdata.v[,2:11] + pdata.v[,16]

### List of data with excess return
data.l <- list(data.e, data.v)

### List of data with real return
pdata.l <- list(pdata.e, pdata.v)

#############################################################################
### Create functions and variables for Chapter 1
#############################################################################

### Create function to generate residual plots
residplotFun <- function(datalist, datalistnum, portfolio){
  data <- datalist[[datalistnum]]
  i <- portfolio + 1
  fit <- lm(data[,i] ~ RMRF + SMB + HML + MOM, data = data)
  
  titlename <- paste(c("(AEW)", "(AVW)")[datalistnum],
                     "Residual Plots: Portfolio", portfolio, sep = " ")
  
  par(mfrow=c(3,2), oma=c(0,0,4,0), mar=c(5.1,4.1,1.1,2.1), las=1, cex=0.7)
  
  plot(data[,"Month"], residuals(fit), 
       main = "Time Series of Residuals",
       xlab = "Month", ylab = "Residuals",
       type = "l", col="blue")
  
  plot(fitted(fit), residuals(fit),
       main = "Fitted vs Residuals",
       xlab = "Fitted", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"RMRF"], residuals(fit),
       main = "RMRF vs Residuals",
       xlab = "RMRF", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"SMB"], residuals(fit),
       main = "SMB vs Residuals",
       xlab = "SMB", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"HML"], residuals(fit),
       main = "HML vs Residuals",
       xlab = "HML", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"MOM"], residuals(fit),
       main = "MOM vs Residuals",
       xlab = "MOM", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  mtext(titlename, 3, 1, outer=TRUE, cex=1)
  
  # Restore default clipping rect
  par(mar = c(5, 4, 4, 2) + 0.1)
  par(mfrow=c(1,1))
}

#############################################################################

### Create function to compute p-value of Lilliefors test for normality
normtestpvalFun <- function(datalist){
  portfolio <- 1:10
  pvalAEW <- c(NA, 10)
  pvalAVW <- c(NA, 10)
  nullAEW <- c(NA, 10)
  nullAVW <- c(NA, 10)
  
  for(i in 1:10){
    fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[1]])
    resid1 <- residuals(fit1)
    pvalAEW[i] <- lillie.test(resid1)$p.value
    if(pvalAEW[i] < 0.05){
      nullAEW[i] <- "Reject"
    } else {nullAEW[i] <- "Accept"}
    
    fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[2]])
    resid2 <- residuals(fit2)
    pvalAVW[i] <- lillie.test(resid2)$p.value
    if(pvalAVW[i] < 0.05){
      nullAVW[i] <- "Reject"
    } else {nullAVW[i] <- "Accept"}
  }
  
  return(data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW))
}

### Compute the p-value
normtestpval.A <- normtestpvalFun(data.l)

#############################################################################

### Create function to compute mean of residuals
meanresidFun <- function(datalist){
  portfolio <- 1:10
  meanresidAEW <- c(NA, 10)
  meanresidAVW <- c(NA, 10)
  
  for(i in 1:10){
    fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[1]])
    meanresidAEW[i] <- mean(residuals(fit1))
    
    fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[2]])
    meanresidAVW[i] <- mean(residuals(fit2))
  }
  
  return(data.frame(portfolio, meanresidAEW, meanresidAVW))
}

### Compute the mean of residuals
meanresid.A <- meanresidFun(data.l)

#############################################################################

### Compute the p-values of non-constant variance score test
datalist <- data.l
portfolio <- 1:10
pvalAEW <- c(NA, 10)
pvalAVW <- c(NA, 10)
nullAEW <- c(NA, 10)
nullAVW <- c(NA, 10)

for(i in 1:10){
  fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[1]])
  pvalAEW[i] <- ncvTest(fit1)$p
  if(pvalAEW[i] < 0.05){
    nullAEW[i] <- "Reject"
  } else {nullAEW[i] <- "Accept"}
  
  fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[2]])
  pvalAVW[i] <- ncvTest(fit2)$p
  if(pvalAVW[i] < 0.05){
    nullAVW[i] <- "Reject"
  } else {nullAVW[i] <- "Accept"}
}

ncvtestpval.A <- data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW)
rm(datalist, portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW, i, fit1, fit2)

#############################################################################
### Create functions and variables for Chapter 2
#############################################################################

### Create function to plot Cook's distance
CiPlotFun <- function(datalist, datalistnum, portfolio, toplot = TRUE){
    data <- datalist[[datalistnum]]
    colnum <- portfolio + 1
    
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    
    ti <- as.matrix(rep(1,N))
    Ci <- as.matrix(rep(1,N))
    for(j in 1:N){
        ti[j] <- resid[j] / (sqrt(sigma2hat) * sqrt(1 - P[j,j]))
        Ci[j] <- ( ((ti[j])^2) / p ) * ( P[j,j] / (1 - P[j,j]) )
    }
    
    calpt.find <- Ci > 1
    Month <- as.matrix(data[,1])
    Cidata <- cbind(Month, Ci)
    colnames(Cidata) <- c("Month", "Ci")
    Cidata <- as.data.frame(Cidata)
    Cidata[,1] <- as.Date(Cidata[,1])
    calpt.count <- sum(calpt.find)
    Cidata.sort <- Cidata[order(-Ci),]
    
    dataname <- c("(AEW)", "(AVW)")[datalistnum]
    if(toplot == TRUE){
        plot(Cidata[,"Month"], Cidata[,"Ci"],
        main = paste(dataname, "\nCook's Distance: Portfolio", portfolio,
        sep = " "),
        xlab = "Month", ylab = "Cook's Distance", col="blue")
        points(Cidata[,"Month"][calpt.find], Cidata[,"Ci"][calpt.find],
        col="red", pch=19)
    } else {
        hist(Cidata[,"Ci"], breaks = 100,
        col = "skyblue", xlab = "Cook's Distance",
        main = paste(dataname, "Histogram\nCook's Distance: Portfolio",
        portfolio, sep = " "))
    }
}

#############################################################################

### Create function to plot covariance ratio
CVRiPlotFun <- function(datalist, datalistnum, portfolio, toplot = TRUE){
    data <- datalist[[datalistnum]]
    colnum <- portfolio + 1
    
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    sigma2hat.i <- as.matrix(rep(1,N))
    for(j in 1:N){
        sigma2hat.i[j] <- (N - p)*sigma2hat/(N - p -1) -
        resid[j]^2/((N - p - 1)*(1 - P[j,j]))
    }
    
    CVRi <- as.matrix(rep(1,N))
    for(j in 1:N){
        CVRi[j] <- ((sigma2hat.i[j]/sigma2hat)^p) / (1 - P[j,j])
    }
    
    calpt.find <- abs(CVRi - 1 ) > 3*p/N
    abs.CVRi_1 <- abs(CVRi - 1 )
    calpt.crit <- 3*p/N
    Month <- as.matrix(data[,1])
    CVRidata <- data.frame(Month, CVRi, abs.CVRi_1, calpt.crit)
    CVRidata[,1] <- as.Date(CVRidata[,1])
    
    calpt.count <- sum(calpt.find)
    CVRidata.sort <- CVRidata[order(-abs.CVRi_1),]
    
    dataname <- c("(AEW)", "(AVW)")[datalistnum]
    if(toplot == TRUE){
        plot(CVRidata$Month, CVRidata$CVRi,
        main = paste(dataname, "\nCovariance Ratio: Portfolio", portfolio,
        sep = " "),
        xlab = "Month", ylab = "Covariance Ratio", col="blue")
        points(CVRidata$Month[calpt.find], CVRidata$CVRi[calpt.find], col="red",
        pch=19)
    } else {
        hist(CVRidata$CVRi, breaks = 100,
        col = "skyblue", xlab = "Covariance Ratio",
        main = paste(dataname, "Histogram\nCovariance Ratio: Portfolio",
        portfolio, sep = " "))
    }
}

#############################################################################

### Create function to combine plots of Cook's distance and covariance ratio
toplotFun <- function(datalist, portfolio, toplot = TRUE){
    if(toplot == TRUE){
        par(mfcol=c(2,2))
        for(i in 1:2){
            CiPlotFun(datalist, i, portfolio)
            CVRiPlotFun(datalist, i, portfolio)
        }
        par(mfrow=c(1,1))
    } else {
        par(mfcol=c(4,1))
        for(i in 1:2){
            CiPlotFun(datalist, i, portfolio, FALSE)
            
        }
        for(i in 1:2){
            CVRiPlotFun(datalist, i, portfolio, FALSE)
        }
        par(mfrow=c(1,1))
    }
}

#############################################################################

### Create function to plot MOM vs monthly portfolio returns
momreturnplotFun <- function(datalist, datalistnum, portfolio){
    data <- datalist[[datalistnum]]
    i <- c(1:7)
    i <- 130 * i
    start <- c(1, i + 1)
    end <- c(i, nrow(data))
    
    j <- portfolio + 1
    temp <- data[,j] + data[,16]
    ymin <- min(min(temp), min(data[,15]))
    ymax <- max(max(temp), max(data[,15]))
    
    titlename <- paste(c("(AEW)", "(AVW)")[datalistnum],
    "MOM & Monthly Portfolio Return: Portfolio", portfolio,
    sep = " ")
    
    par(mfrow=c(4,2), oma=c(3,3,4,0), mar=c(4,2,1,1), las=1, cex=0.7)
    for(k in 1:8){
        plot(data[start[k]:end[k],1], temp[start[k]:end[k]],
        xlab = "", ylab = "", ylim=c(ymin,ymax), type="l", col = "blue")
        lines(data[start[k]:end[k],1], data[start[k]:end[k],15], type="l",
        col = "red")
        legend("topleft", c("Return", "MOM"), cex=0.5, col=c("blue", "red"),
        lty=1, bty="n")
    }
    mtext("Month", 1, 1, outer=TRUE, cex = 0.75)
    mtext("MOM & Monthly Portfolio Return", 2, 1, outer=TRUE, las=0,
    cex = 0.75)
    mtext(titlename, 3, 1, outer=TRUE, cex=1)
    
    # Restore default clipping rect
    par(mar = c(5, 4, 4, 2) + 0.1)
    par(mfrow=c(1,1))
}

#############################################################################

### Create function to compute the Cook's distance for Influential Points
CiCalptFun <- function(data, colnum){
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    
    ti <- as.matrix(rep(1,N))
    Ci <- as.matrix(rep(1,N))
    for(j in 1:N){
        ti[j] <- resid[j] / (sqrt(sigma2hat) * sqrt(1 - P[j,j]))
        Ci[j] <- ( ((ti[j])^2) / p ) * ( P[j,j] / (1 - P[j,j]) )
    }
    
    calpt.find <- Ci > 1
    Month <- as.matrix(data[,1])
    Cidata <- cbind(Month, Ci)
    colnames(Cidata) <- c("Month", "Ci")
    Cidata <- as.data.frame(Cidata)
    Cidata[,1] <- as.Date(Cidata[,1])
    calpt.count <- sum(calpt.find)
    Cidata.sort <- Cidata[order(-Ci),]
    
    if(calpt.count == 0){
        return(NULL)
    } else {
        return(Cidata.sort[1:calpt.count,])
    }
    
}

#############################################################################

### Create function to compute the covariance ratios for Influential Points
CVRiCalptFun <- function(data, colnum){
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    sigma2hat.i <- as.matrix(rep(1,N))
    
    for(j in 1:N){
        sigma2hat.i[j] <- (N - p)*sigma2hat/(N - p -1) -
        resid[j]^2/((N - p - 1)*(1 - P[j,j]))
    }
    
    CVRi <- as.matrix(rep(1,N))
    for(j in 1:N){
        CVRi[j] <- ((sigma2hat.i[j]/sigma2hat)^p) / (1 - P[j,j])
    }
    
    calpt.find <- abs(CVRi - 1 ) > 3*p/N
    abs.CVRi_1 <- abs(CVRi - 1 )
    calpt.crit <- 3*p/N
    Month <- as.matrix(data[,1])
    CVRidata <- data.frame(Month, CVRi, abs.CVRi_1, calpt.crit)
    CVRidata[,1] <- as.Date(CVRidata[,1])
    
    calpt.count <- sum(calpt.find)
    CVRidata.sort <- CVRidata[order(-abs.CVRi_1),]
    
    if(calpt.count == 0){
        return(NULL)
    } else {
        return(CVRidata.sort[1:calpt.count,])
    }
}

#############################################################################

### Create function to compute influence measures of influential observations
exclusioninfoFun <- function(data, portfolio){
    Ci <- as.list(c(1:10))
    CVARi <- as.list(c(1:10))
    
    for(i in 1:10){
        
        if(is.null(CiCalptFun(data, i+1)) == TRUE){
            Ci[[i]] <- "None"
        } else {
            Ci[[i]] <- CiCalptFun(data, i+1)
        }
        
        if(is.null(CVRiCalptFun(data, i+1)) == TRUE){
            CVARi[[i]] <- "None"
        } else {
            CVARi[[i]] <- CVRiCalptFun(data, i+1)
        }
        
    }
    
    exclusion.info <- list("Cook's Distance" = Ci[[portfolio]],
    "Covariance Ratio" = CVARi[[portfolio]])
    
    return(exclusion.info)
}

#############################################################################
### Create functions and variables for Chapter 3
#############################################################################

### Get the excluded months
exclusion.Be <- c()
exclusion.Bv <- c()
exclusion.Ce <- c()
exclusion.Cv <- c()

for(i in 1:10){
    exclusion.Be <- c(exclusion.Be,
    CiCalptFun(data.e, i+1)[,1])
    exclusion.Bv <- c(exclusion.Bv,
    CiCalptFun(data.v, i+1)[,1])
    exclusion.Ce <- c(exclusion.Ce,
    CVRiCalptFun(data.e, i+1)[,1])
    exclusion.Cv <- c(exclusion.Cv,
    CVRiCalptFun(data.v, i+1)[,1])
}

exclusion.Be <- as.Date(sort(unique(exclusion.Be)))
exclusion.Bv <- as.Date(sort(unique(exclusion.Bv)))
exclusion.Ce <- as.Date(sort(unique(exclusion.Ce)))
exclusion.Cv <- as.Date(sort(unique(exclusion.Cv)))

### Get the subset lists
subset.Ae <- as.list(c(1:10))
subset.Av <- as.list(c(1:10))
subset.Be <- as.list(c(1:10))
subset.Bv <- as.list(c(1:10))
subset.Ce <- as.list(c(1:10))
subset.Cv <- as.list(c(1:10))

for(i in 1:10){
    subset.Ae[[i]] <- rep(TRUE, nrow(data.e))
    subset.Av[[i]] <- rep(TRUE, nrow(data.v))
    subset.Be[[i]] <- !(data.e[,1] %in% exclusion.Be)
    subset.Bv[[i]] <- !(data.v[,1] %in% exclusion.Bv)
    subset.Ce[[i]] <- !(data.e[,1] %in% exclusion.Ce)
    subset.Cv[[i]] <- !(data.v[,1] %in% exclusion.Cv)
}

### Panel B data
  data.Be <- data.e[subset.Be[[1]],] ### AEW excess return
  data.Bv <- data.v[subset.Bv[[1]],] ### AVW excess return
  pdata.Be <- pdata.e[subset.Be[[1]],] ### AEW real return
  pdata.Bv <- pdata.v[subset.Bv[[1]],] ### AVW real return
  
  ### List of data
  data.Bl <- list(data.Be, data.Bv)
  pdata.Bl <- list(pdata.Be, pdata.Bv)

### Panel C data 
  data.Ce <- data.e[subset.Ce[[1]],] ### AEW excess return
  data.Cv <- data.v[subset.Cv[[1]],] ### AVW excess return
  pdata.Ce <- pdata.e[subset.Ce[[1]],] ### AEW real return
  pdata.Cv <- pdata.v[subset.Cv[[1]],] ### AVW real return
  
  ### List of data
  data.Cl <- list(data.Ce, data.Cv)
  pdata.Cl <- list(pdata.Ce, pdata.Cv)

#############################################################################

### Create function to compute the statistics
olsFun <- function(data, portfolio){
  varindex <- portfolio + 1
  fit <- lm(data[,varindex] ~ RMRF + SMB + HML + MOM, data=data)
  
  varnames <- c("Intercept", "RMRF", "SMB", "HML", "MOM")
  beta <- matrix(NA, 1, 5)
  beta[1,] <- fit$coefficients
  colnames(beta) <- varnames
  
  beta.t <- matrix(NA, 1, 5)
  beta.t[1,] <- coeftest(fit)[,"t value"]
  colnames(beta.t) <- varnames
  
  beta.NWt <- matrix(NA, 1, 5)
  beta.NWt[1,] <- coeftest(fit, NeweyWest(fit))[,"t value"]
  colnames(beta.NWt) <- varnames
  
  return(list("Beta" = beta, "Beta t-statistic" = beta.t,
              "Beta Newey-West t-statistic" = beta.NWt))
}

#############################################################################

### Create function to generate the panel summary
panelsummaryFun <- function(data, subset){
    panel.sum <- olsFun(data[subset[[1]], ], 1)
    for(i in 2:10){
        stats.temp <- olsFun(data[subset[[i]], ], i)
        for(j in 1:3){
            panel.sum[[j]] <- rbind(panel.sum[[j]], stats.temp[[j]])
        }
    }
    
    Portfolio <- 1:10
    for(j in 1:3){
        panel.sum[[j]] <- cbind(Portfolio, panel.sum[[j]])
    }
    
    return(panel.sum)
}

#############################################################################

### Generate summary for each panel
### AEW
panelAe.summary <- panelsummaryFun(data.e, subset.Ae)
panelBe.summary <- panelsummaryFun(data.e, subset.Be)
panelCe.summary <- panelsummaryFun(data.e, subset.Ce)

### AVW
panelAv.summary <- panelsummaryFun(data.v, subset.Av)
panelBv.summary <- panelsummaryFun(data.v, subset.Bv)
panelCv.summary <- panelsummaryFun(data.v, subset.Cv)

#############################################################################

### (AEW) Beta
### Panel A
beta.Ae <- panelAe.summary[[1]]

### Panel B
beta.Be <- panelBe.summary[[1]]

### Percentage change from panel A to panel B
beta.ABe <- ((panelBe.summary[[1]] - panelAe.summary[[1]]) /
panelAe.summary[[1]]) * 100
beta.ABe[,1] <- 1:10

#############################################################################

### (AVW) Beta
### Panel A
beta.Av <- panelAv.summary[[1]]

### Panel B
beta.Bv <- panelBv.summary[[1]]

### Percentage change from panel A to panel B
beta.ABv <- ((panelBv.summary[[1]] - panelAv.summary[[1]]) /
panelAv.summary[[1]]) * 100
beta.ABv[,1] <- 1:10

#############################################################################

### (AEW) Beta t-statistic
### Panel A
beta_t.Ae <- panelAe.summary[[2]]

### Panel C
beta_t.Ce <- panelCe.summary[[2]]

### Percentage change from panel A to panel C
beta_t.ACe <- ((panelCe.summary[[2]] - panelAe.summary[[2]]) /
panelAe.summary[[2]]) * 100
beta_t.ACe[,1] <- 1:10

#############################################################################

### (AEW) Beta Newey-West t-statistic
### Panel A
beta_nwt.Ae <- panelAe.summary[[3]]

### Panel C
beta_nwt.Ce <- panelCe.summary[[3]]

### Percentage change from panel A to panel C
beta_nwt.ACe <- ((panelCe.summary[[3]] - panelAe.summary[[3]]) /
panelAe.summary[[3]]) * 100
beta_nwt.ACe[,1] <- 1:10

#############################################################################

### (AVW) Beta t-statistic
### Panel A
beta_t.Av <- panelAv.summary[[2]]

### Panel C
beta_t.Cv <- panelCv.summary[[2]]

### Percentage change from panel A to panel C
beta_t.ACv <- ((panelCv.summary[[2]] - panelAv.summary[[2]]) /
panelAv.summary[[2]]) * 100
beta_t.ACv[,1] <- 1:10

#############################################################################

### (AVW) Beta Newey-West t-statistic
### Panel A
beta_nwt.Av <- panelAv.summary[[3]]

### Panel C
beta_nwt.Cv <- panelCv.summary[[3]]

### Percentage change from panel A to panel C
beta_nwt.ACv <- ((panelCv.summary[[3]] - panelAv.summary[[3]]) /
panelAv.summary[[3]]) * 100
beta_nwt.ACv[,1] <- 1:10

#############################################################################

### Examine the OLS assumptions for panels B and C

### Panel B

### Compute the p-Values of Lilliefors test for normality
normtestpval.B <- normtestpvalFun(data.Bl)

### Compute the p-values of non-constant variance score test
datalist <- data.Bl
portfolio <- 1:10
pvalAEW <- c(NA, 10)
pvalAVW <- c(NA, 10)
nullAEW <- c(NA, 10)
nullAVW <- c(NA, 10)

for(i in 1:10){
  fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[1]])
  pvalAEW[i] <- ncvTest(fit1)$p
  if(pvalAEW[i] < 0.05){
    nullAEW[i] <- "Reject"
  } else {nullAEW[i] <- "Accept"}
  
  fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[2]])
  pvalAVW[i] <- ncvTest(fit2)$p
  if(pvalAVW[i] < 0.05){
    nullAVW[i] <- "Reject"
  } else {nullAVW[i] <- "Accept"}
}

ncvtestpval.B <- data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW)
rm(datalist, portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW, i, fit1, fit2)

### Panel C

### Compute the p-Values of Lilliefors test for normality
normtestpval.C <-normtestpvalFun(data.Cl)

### Compute the p-values of non-constant variance score test
datalist <- data.Cl
portfolio <- 1:10
pvalAEW <- c(NA, 10)
pvalAVW <- c(NA, 10)
nullAEW <- c(NA, 10)
nullAVW <- c(NA, 10)

for(i in 1:10){
  fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[1]])
  pvalAEW[i] <- ncvTest(fit1)$p
  if(pvalAEW[i] < 0.05){
    nullAEW[i] <- "Reject"
  } else {nullAEW[i] <- "Accept"}
  
  fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[2]])
  pvalAVW[i] <- ncvTest(fit2)$p
  if(pvalAVW[i] < 0.05){
    nullAVW[i] <- "Reject"
  } else {nullAVW[i] <- "Accept"}
}

ncvtestpval.C <- data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW)
rm(datalist, portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW, i, fit1, fit2)


#############################################################################
### Create functions and variables for Chapter 4
#############################################################################

### Create function to compute mean and sd of monthly portfolio return

statFun <- function(pdata, subset){
    statsummary <- matrix(NA, 10, 3)
    rownames(statsummary) <- rep("", 10)
    colnames(statsummary) <- c("Portfolio", "Mean", "SD")
    statsummary[,1] <- 1:10
    for(i in 1:10){
        statsummary[i,2] <- mean(pdata[subset[[i]], i+1])
        statsummary[i,3] <- sd(pdata[subset[[i]], i+1])
    }
    return(statsummary)
}

#############################################################################

### (AEW) Mean and sd of return
### Panel A
meansd.Ae <- statFun(pdata.e, subset.Ae)

### Panel B
meansd.Be <- statFun(pdata.e, subset.Be)

### Panel C
meansd.Ce <- statFun(pdata.e, subset.Ce)

### Percentage change from panel A to panel B
meansd.ABe <- ((meansd.Be - meansd.Ae) / meansd.Ae) * 100
meansd.ABe[,1] <- 1:10

### Percentage change from panel A to panel C
meansd.ACe <- ((meansd.Ce - meansd.Ae) / meansd.Ae) * 100
meansd.ACe[,1] <- 1:10

#############################################################################

### (AVW) Mean and sd of return
### Panel A
meansd.Av <- statFun(pdata.v, subset.Av)

### Panel B
meansd.Bv <- statFun(pdata.v, subset.Bv)

### Panel C
meansd.Cv <- statFun(pdata.v, subset.Cv)

### Percentage change from panel A to panel B
meansd.ABv <- ((meansd.Bv - meansd.Av) / meansd.Av) * 100
meansd.ABv[,1] <- 1:10

### Percentage change from panel A to panel C
meansd.ACv <- ((meansd.Cv - meansd.Av) / meansd.Av) * 100
meansd.ACv[,1] <- 1:10

#############################################################################

### Create function to compute numbers for mean-variance optimal portfolio

effpfFun <- function(pdata, mu.min, mu.max, nobs, toshorts = TRUE){
    
    pfoptFun <- function(mu, pdata, toshorts){
        pfmu <- matrix(colMeans(pdata[,2:11]),nrow=1)
        pfcov <- cov(pdata[,2:11])
        if(toshorts == TRUE){
            pfopt.sol <- portfolio.optim(x = pfmu, pm = mu,
            covmat = pfcov, shorts = toshorts)
        } else {
            pfopt.sol <- portfolio.optim(x = pfmu, pm = mu,
            covmat = pfcov, shorts = toshorts,
            reslow = rep(0,10), reshigh = rep(1,10))
        }
        
        wt <- pfopt.sol$pw
        SD <- sqrt((t(wt) %*% pfcov %*% wt)[1])
        return(c(mu, SD, wt))
    }
    
    mu.list <- seq(mu.min, mu.max, length.out = nobs)
    effpf <- t(sapply(mu.list, pfoptFun, pdata, toshorts))
    colnames(effpf) <- c("Mean", "SD", paste("W", 1:10, sep=""))
    rownames(effpf) <- rep("", nrow(effpf))
    
    return(effpf)
}

#############################################################################

### (AEW) Mean-Variance Optimal Portfolio

### Mean of each monthly portfolio return
pfmu.e <- matrix(colMeans(pdata.e[,2:11]),nrow=1)

### Max and min of the mean
pmin.e <- min(pfmu.e) + 1e-4
pmax.e <- max(pfmu.e) - 1e-4

### Mean-Variance Optimal Portfolio
effpfT.e <- effpfFun(pdata.e, 0, 0.5, 1000, TRUE)
effpfF.e <- effpfFun(pdata.e, pmin.e, pmax.e, 100, FALSE)

### Subset for nicer layout
min_x.e <- min(effpfT.e[,2], effpfF.e[,2])
max_x.e <- max(effpfF.e[,2])
effpfT.e <- effpfT.e[effpfT.e[,2] <= max_x.e,]

#############################################################################

### (AVW) Mean-Variance Optimal Portfolio

### Mean of each monthly portfolio return
pfmu.v <- matrix(colMeans(pdata.v[,2:11]),nrow=1)

### Max and min of the mean
pmin.v <- min(pfmu.v) + 1e-4
pmax.v <- max(pfmu.v) - 1e-4

### Mean-Variance Optimal Portfolio
effpfT.v <- effpfFun(pdata.v, 0, 0.5, 1000, TRUE)
effpfF.v <- effpfFun(pdata.v, pmin.v, pmax.v, 100, FALSE)

### Subset for nicer layout
min_x.v <- min(effpfT.v[,2], effpfF.v[,2])
max_x.v <- max(effpfF.v[,2])
effpfT.v <- effpfT.v[effpfT.v[,2] <= max_x.v,]

#############################################################################

### (Shorts Allowed) Mean-Variance Optimal Portfolio: AEW vs AVW
pmin <- 0
pmax <- 0.03

effpf.e <- effpfFun(pdata.e, pmin, pmax, 31, TRUE)
effpf.v <- effpfFun(pdata.v, pmin, pmax, 31, TRUE)

min_x.ev <- min(effpf.e[,2], effpf.v[,2])
max_x.ev <- max(effpf.e[,2], effpf.v[,2])

#############################################################################
### The End of R codes for functions and variables
#############################################################################
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract}

% The text of the abstract goes here.  If you need to use a \section
% command you will need to use \section*, \subsection*, etc. so that
% you don't get any numbering.  You probably won't be using any of
% these commands in the abstract anyway.

We examine whether OLS is the right method to apply on the Fama French Carhart Four Model using the observed data of the four-factor and the monthly returns of each of the 10 portfolios formed monthly on momentum, where the 10 portfolios are constructed monthly using NYSE prior (2-12) return in ascending order with decile breakpoints. Next, we use Cook's distance and covariance ratio to measure the influence of an observation on betas and t-statistics respectively. Finally, we compare two mean-variance optimal portfolios, in which the portfolios are both comprised of the 10 portfolios formed on momentum but with two different weighted methods.

\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Introduction}

\section{The Fama French Carhart Four Model}
The Fama French Carhart Four Model \cite{mom1} is a 4-factor model using Fama and French's 3-factor model plus an additional factor capturing Jegadeesh and Titman's one-year momentum anomaly. The Fama French Carhart Four Model is formulated as
\begin{equation} 
\label{eq:eq1.1}
r_{i,t} = \alpha_i + \beta_{i,RMRF}RMRF_{t} + \beta_{i,SMB}SMB_t + \beta_{i,HML}HML_t + \beta_{i,MOM}MOM_t + \epsilon_{i,t}
\end{equation}
where i = 1,...,10 indexes the 10 portfolios formed monthly on momentum factor, t indexes for month, r is the portfolio return in excess of the risk-free rate, RMRF, SMB, and HML are the Fama-French market, size, and value factors, MOM is the momentum factor, $\epsilon$ is the unobservable random error, and the various $\beta$s are the corresponding factor loadings. The paper focus on the empirical study of the equation \eqref{eq:eq1.1} using the data from Kenneth R. French's Data Library. \footnote{The Kenneth R. French's Data Library \url{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html} contains data for U.S. research returns, historical benchmark returns, U.S. research breakpoints, U.S. book equity, international research returns, and developed market factors and returns.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Data}
We are using the following three sets of data from Kenneth R. French's Data Library \cite{data1}.
\begin{itemize}
  \item Fama/French Factors
  \item Momentum Factor (MOM)
  \item 10 Portfolios Formed Monthly on Momentum
\end{itemize}

\subsection{Fama/French Factors}
This dataset includes the monthly data for the three Fama/French Factors RMRF, SMB, and HML. It also includes the risk-free rate factor $R_f$ so that we can get the market rate factor $R_m$ by adding the risk-free factor to RMRF factor.
\\

SMB (Small Minus Big) \cite{data1} is the average return on the three small portfolios minus the average return on the three big portfolios, which is equal to
\begin{center}
$\frac{1}{3}(Small \, Value + Small \, Neutral + Small \, Growth) - \frac{1}{3}(Big \, Value + Big \, Neutral + Big \, Growth)$
\end{center}
\\

HML (High Minus Low) \cite{data1} is the average return on the two value portfolios minus the average return on the two growth portfolios, which is equal to
\begin{center}
$\frac{1}{2}(Small \, Value + Big \, Value) - \frac{1}{2}(Small \, Growth + Big \, Growth)$
\end{center}
\\

RMRF ($R_m - R_f$) \cite{data1} is the excess return on the market, value-weight return of all CRSP firms incorporated in the US and listed on the NYSE, AMEX, or NASDAQ that have a CRSP share code of 10 or 11 at the beginning of month t, good shares and price data at the beginning of t, and good return data for t minus the one-month Treasury bill rate.
\\

These Fama/French factors \cite{data1} are constructed using the 6 value-weight portfolios formed on size and book-to-market. The portfolios, which are constructed at the end of each June, are the intersections of 2 portfolios formed on size (market equity, ME) and 3 portfolios formed on the ratio of book equity to market equity (BE/ME). The size breakpoint for year t is the median NYSE market equity at the end of June of year t. BE/ME for June of year t is the book equity for the last fiscal year end in $t - 1$ divided by ME for December of $t - 1$. The BE/ME breakpoints are the 30th and 70th NYSE percentiles.

\subsection{Momentum Factor (MOM)}
This dataset includes the monthly data for the Momentum (MOM) factor \cite{data1}, MOM is the average return on the two high prior return portfolios minus the average return on the two low prior return portfolios, which is equal to
\begin{center}
$\frac{1}{2}(Small \, High + Big \, High) - \frac{1}{2}(Small \, Low + Big \, Low)$
\end{center}

\subsection{10 Portfolios Formed Monthly on Momentum}
This dataset includes the monthly data for return of 10 portfolios formed monthly on momentum.
\\

The portfolios \cite{data1} constructed each month include NYSE, AMEX, and NASDAQ stocks with prior return data. To be included in a portfolio for month t (formed at the end of month $t - 1$), a stock must have a price for the end of month $t - 13$ and a good return for $t - 2$.
\\

Stocks are sorted at the beginning of the first trading day of each month from January 1925 to August 2013 into decile portfolios based on their previous calendar year's return. We have two sets of data; one is average equal weighted (AEW) and the other one is average value weighted (AVW). The portfolios are weighted monthly so the weights are readjusted whenever a stock disappears. Stocks with the lowest NYSE prior (2-12) return comprise decile 1 and stocks with the highest comprise decile 10 \cite{data1}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Programming Language and Reproducibility}
We use the open-source programming language R for the computation and data analysis that are necessary for our research. All the figures and other output appeared in this thesis are reproducible using the R codes in appendix G together with the three sets of data from Kenneth R. French's Data Library. We will also use GitHub to facilitate reproducible computational results and share our research. All our datasets and programming files, together with this thesis, can be found here:
\url{https://github.com/henrykmwong}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Multiple Linear Regression}
In this chapter, we set up a multiple linear regression model for the equation \eqref{eq:eq1.1} and identify the influential observations for each of the portfolios. The model \cite{model1} is:
  \begin{equation} 
\label{eq:eq1.2}
Y = X\beta + \epsilon
\end{equation}
On the left, $Y$ is an $N \times 1$ vector of observable random variables. The $Y$ vector is called the \textit{dependent} or \textit{response} variable. On the right, the first term is $X$, an $N \times p$ matrix of observable random variables, one column for each data variable. The $X$ matrix is called the \textit{design matrix} and the data variables are called the \textit{explanatory} or \textit{independent} variables. We assume that $N > p$ and the design matrix has full rank. The second term on the right is $\beta$, a $p \times 1$ vector of unknown coefficents or parameters to be estimated. The final term on the right is $\epsilon$, an $N \times 1$ vector of unobservable random errors.
\\

In the application of the Fama French Carhart Four Model on our data, $N$ is the number of monthly observations, which is equal to 1,039 and p is the number of coefficient variables, which is equal to 5 (intercept plus coefficients for the four factors).
\\

We will write $Y_i$ for the $i$th observation of $Y$, $X_i$ for the $i$th row of $X$, and $\epsilon_i$ for the ith component of $\epsilon$. The matrix equation \eqref{eq:eq1.2} is equivalent to a set of $N$ ordinary equations, one for each observation. For the $i$th observation, the equation is 
\begin{equation} 
\label{eq:eq1.3}
Y_i = X_{i}\beta + \epsilon_i
\end{equation}

To estimate the $\beta$ in equation \eqref{eq:eq1.2}, we use the method of ordinary least squares (OLS) with the following assumptions \cite{model1}.
\begin{enumerate}
\item The $\epsilon_i$ are independent and identically distributed with mean zero and constant variance (or homoscedasticity) $\sigma^2$.
\item If $X$ is random, $\epsilon$ is independent of $X$.
\end{enumerate}
\\

Using the method of OLS, we have:
\begin{equation} 
\label{eq:eq1.4}
\hat{\beta} = (X^{T}X)^{-1}X^{T}Y
\end{equation}

\begin{equation} 
\label{eq:eq1.5}
\hat{Y} = X\hat{\beta}
\end{equation}

where

\begin{equation} 
\label{eq:eq1.6}
P = X(X^{T}X)^{-1}X^{T}
\end{equation}

\begin{equation} 
\label{eq:eq1.7}
e = Y - \hat{Y} = (I - P)Y
\end{equation}

\begin{equation} 
\label{eq:eq1.8}
{\hat{\sigma}}^{2} = \frac{e^{T}e}{N - p}
\end{equation}

\begin{equation} 
\label{eq:eq1.9}
\begin{align*}
\Psi &= cov(\hat{\beta}|X)\\
&= \sigma^{2}(X^{T}X)^{-1}
\end{align*}
\end{equation}

\begin{equation} 
\label{eq:eq1.10}
cov(\hat{Y}|X) = \sigma^{2}P
\end{equation}

and

\begin{equation} 
\label{eq:eq1.11}
cov(e|X) = \sigma^{2}(I - P)
\end{equation}

We denote the vector of estimeated coefficients and the vector of fitted values by $\hat{\beta}$ and $\hat{Y}$. The $N \times N$ matrix $P$ in equation \eqref{eq:eq1.6} is called the \textit{hat matrix}. It is also called the \textit{projection matrix} because it projects $Y$ into the column space of $X$. The vector of the (least squares) residuals is denoted by $e$ in equation \eqref{eq:eq1.7}.
\\

The following equation will be used in this chapter for identifying influential observations.

\begin{equation} 
\label{eq:eq1.12}
t_i = \frac{e_i}{\hat{\sigma}\sqrt{1 - p_i}}
\end{equation}

where $e_i$ is the $i$th component of the vector of residuals $e$ and $p_i$ is the $i$th diagonal element of P.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Regression Diagnostic}
The OLS assumptions are difficult to check because $\epsilon$ is not observable \cite{model1}. One way to work around the difficulty is to use the OLS residuals. 
\\

We can use the Lilliefors test \cite{norm1} to test the null hypothesis that the residuals come from the Gaussian distribution. The 2nd and 3rd columns of table \ref{table:t1.1} show the p-value of the Lilliefors test for each portfolio from each weighting method. We reject all the null that the residuals are Gaussian distributed with $5\%$ $\alpha$ level and therefore the OLS estimates are not equivalent to the maximum-likelihood estimates. However, it does not violate the first OLS assumption as long as the residuals are independent and identically distributed with mean zero and constant variance.
\\

<<echo=FALSE, results="asis">>=
print(
  xtable(normtestpval.A, 
         caption = "Lilliefors Test for Normality with 0.05 Alpha Level",
         label = "table:t1.1", display = c(rep('d', 2), rep('e', 2), rep('s', 2))),
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

As shown in appendix A, each of the scatterplots of residuals versus the observed values of the independent variables does not show any pattern, therefore the second OLS assumption seem to be valid. For illustration, figure \ref{fig:f1.1} shows the residual plots for portfolio 2 with average equal weighted (AEW) method. In addition, as shown in table \ref{table:t1.2}, the mean of residuals for each portfolio and weighting method is equal to zero with accuracy of 18 decimal places.  
\\

<<echo=FALSE, results="asis">>=
print(
  xtable(meanresid.A, 
         caption = "Mean of Residuals",
         label = "table:t1.2", display = c(rep('d', 2), rep('e', 2))), 
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

However, based on the times series of the residuals in appendix A, the residuals seem not to be independent and identically distributed because the residuals have larger fluctuations before 1940 and after 1980. The time series plots also suggest that the homoscedasticity (constant variance) assumption is violated. We can use the non-constant variance score test \cite{var1} to test the null hypothesis that the residuals are distributed with constant variance. The 2nd and 3rd columns of table \ref{table:t1.3} shows the p-value of the non-constant variance score test for each portfolio and weighting method. We reject all the null that the residuals are distributed with constant variance with $5\%$ $\alpha$ level except for portfolios 7 and 10 with average value weighted method.
\\

Our failure of meeting some of the OLS assumptions suggests that our dataset constains influential observations. Chapter 2 will discuss the identification of the influential observations. Even though we do not meet all the OLS assumptions, we can still apply the method of OLS on the multiple linear regression model for the equation \eqref{eq:eq1.1}, but we need to exercise care in interpreting the results. In particular, we can use Newey-West t-statistic to handle the autocorrelation and heteroscedasticity problem, which will be discussed next.
\\

<<echo=FALSE, results="asis">>=
print(
  xtable(ncvtestpval.A, 
         caption = "Non-constant Variance Score Test with 0.05 Alpha Level",
         label = "table:t1.3", display = c(rep('d', 2), rep('e', 2), rep('s', 2))),
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

\begin{figure}[h!]
\caption{Residual Plots (e.g. AEW Portfolio 2)}
\label{fig:f1.1}
<<echo=FALSE>>=
residplotFun(data.l, 1, 2)
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Heteroskedasticity and Autocorrelation Consistent Covariance}
When random errors exhibit autocorrelation and/or heteroskedasticity, $cov(\epsilon|X)$ will not be equal to $\sigma^{2}I$, therefore equation \eqref{eq:eq1.9} no longer holds. In other words, $\hat{\sigma}^{2}(X^{T}X)^{-1}$, denoted by $\hat{\Psi}$, will be a biased estimator for $cov(\hat{\beta}|X)$. Instead, we have\\

\begin{equation} 
\label{eq:eq1.13}
E(\epsilon|X) = 0_{N \times 1}
\end{equation}

\begin{equation} 
\label{eq:eq1.14}
cov(\epsilon|X) = \Omega
\end{equation}

where $\Omega$ is a $N \times N$ matrix. The formula for $cov(\hat{\beta}|X)$ becomes:\\

\begin{equation} 
\label{eq:eq1.15}
\begin{align*}
cov(\hat{\beta}|X) &= (X^{T}X)^{-1}X^{T} \Omega X(X^{T}X)^{-1}\\
&= (\frac{1}{N}X^{T}X)^{-1} \frac{1}{N}\Phi (\frac{1}{N}X^{T}X)^{-1}
\end{align*}
\end{equation}

where $\Phi = N^{-1} X^{T} \Omega X$.\\

To tackle the biased problem, Newey and West \cite{tstat1} have suggested to use the heteroskedasticity and autocorrelation consistent (HAC) covariance estimator $\hat{\Psi}_{HAC}$, which can be computed by substituting $\hat{\Phi}$ into equation \eqref{eq:eq1.15} with 

\begin{equation} 
\label{eq:eq1.16}
\hat{\Phi} = \frac{1}{N}\sum_{i,j=1}^{N} \omega_{|i-j|} \hat{V}_i {\hat{V}_j}^T
\end{equation}

where $\omega = (\omega_0,\ldots,\omega_{N-1})^T$ is a vector of weight with lag $l = |i-j|$. Furthermore, Newey and West \cite{tstat1} have proved that $\hat{\Psi}_{HAC}$ is positive semi-definite. Thus, in the presence of heteroskedasticity and autocorrelation, to test the null hypothesis that $\beta_k = 0$ against the alternative $\beta_k \neq 0$ with $k = 1,\ldots,p$, we should look at the Newey-West t-statistic:\\

\begin{equation} 
\label{eq:eq1.17}
t_{NW} = \frac{\hat{\beta}_k}{\hat{SE}_{NW}(\beta_k)}
\end{equation}

where $\hat{SE}_{NW}(\beta_k)$ is the $kk$th element of $\hat{\Psi}_{HAC}$, in addition to the ordinary t-statistic.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Identification of Influential Observation}

\section{Influence Measure}

Belsley, Kuh, and Welsch \cite{influ2} have given an appropriate definition of influential observation.

\begin{quote}
An influential observation is one which, either individually or together with several other observations, has a demonstrably larger impact on the calculated values of various estimates (coefficients, standard errors, t-values, etc.) than is the case for most of other observations. One obvious means for examining such an impact is to delete each row, one at a time, and note the resultant effect on the various calculated values. Rows whose deletion produces relatively large changes in the calculated values are deemed influential.
\end{quote}

There are numerous influence measures and there is a great deal of redundancy among the measures \cite{influ1}. After experiencing with several data sets, Chatterjee and Hadi \cite{influ1} have concluded that only a few influence measures are sufficient for identifying influential observations. Influence measures can be classified into the following five groups \cite{influ1}:
\begin{enumerate}
  \item Measures based on residuals
  \item Measures based on the prediction matrix
  \item Measures based on the volume of confidence ellipsoids
  \item Measures based on influence function
  \item Measures based on on partial influence
\end{enumerate}
\\

In particular, measures based on the volume of confidence ellipsoids measure the impact of influential observations on the t-statistic of regression coefficient; whereas the measures based on influence function measure the impact of influential observations on the estimated regression coefficient. For instance, Anderson, Bianchi, and Goldberg \cite{advisor1} use covariance ratio to examine the impact of the influential observations on the Newey-West t-statistic of the FVIX betas. In this thesis, we will use \textit{Cook's Distance} and \textit{Covariance Ratio}, in which Cook's Distance is the measure based on influence function and Covariance Ratio is the measure based on the volume of confidence ellipsoid.

\subsection{Cook's Distance}
Cook's disance was proposed by Cook \cite{cook1} in 1977.
\begin{equation} 
\label{eq:eq2.1}
C_i = \frac{t_i}{p}\frac{p_i}{1 - p_i}
\end{equation}

where $t_i$ is defined in equation \eqref{eq:eq1.12} and $p_i$ is the $i$th diagonal element of P as defined in equation \eqref{eq:eq1.6}.\\

We denote $C_i$ as the Cook's distance for the $i$th observation. Chatterjee and Hadi \cite{influ1} interpret $C_i$ as:
\begin{quote}
The scaled Euclidean distance between the two vectors of fitted values when the fitting is done by including or excluding the $i$th observation.
\end{quote}

Cook \cite{cook2} suggests that the $i$th observation can be considered as influential if $C_i$ is greater than 1.

\subsection{Covariance Ratio}
Covariance Ratio was proposed by Belsley, Kuh, and Welsch \cite{influ2} in 1980.
\begin{equation} 
\label{eq:eq2.2}
CVR_i = \left(\frac{N - p - {t_i}^2}{N - p - 1} \right)^p \middle/ (1 - p_i)
\end{equation}
We denote $CVR_i$ as the covariance ratio for the $i$th observation. The data point with $|CVR_i - 1| > 3p/N$ is considered as influential observation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Influential Observation}

With the equations \eqref{eq:eq2.1} and \eqref{eq:eq2.2}, along with the influential point criteria, we can easily find out the influential points for each of the ten portfolios and weighting methods based on the influence measures of Cook's distance and covariance ratio. Appendix B includes the plots of Cook's distance and covariance ratios for each of the ten portfolios and weighting methods. The influential observations for the corresponding influence measures are highlighted with solid red. Appendix D shows the information of the influential observations identified by the two influence measures for each portfolio and weighting methods. 
\\

The plots in appendix B show that that covariance ratio identifies more influential observations than Cook's distance. Figure \ref{fig:f2.1} shows an example for portfolio 2. We then combine our identification results from the ten portfolios. Out of 1,039 total monthly observations, we find that for the average equal weighted method, Cook's distance and covariance ratio identify 5 and 187 influential observations respectively. For the average value weighted method, Cook's distance and covariance ratio identify 4 and 232 influential observations respectively. Appendix E shows the corresponding months of the influential observations.
\\

\begin{figure}[h!] 
\caption{Plot of Cook's Distance and Covariance Ratio (e.g. Portfolio 2)}
\label{fig:f2.1}
<<echo=FALSE>>=
toplotFun(data.l, 2, TRUE)
@
\end{figure}

\begin{figure}[h!] 
\caption{Histogram of Cook's Distance and Covariance Ratio (e.g. Portfolio 2)}
\label{fig:f2.2}
<<echo=FALSE>>=
toplotFun(data.l, 2, FALSE)
@
\end{figure}

\subsection{Influential Observations Identified by Cook's Distance}

Based on the combined identification results in appendix E, the average equal weighted method has one additional influential observation \textbf{1939-09} and the rest of four influential observations are the same as the ones for the average value weighted method. The influential observations identified by Cook's distance are as follows. 
\begin{itemize}
  \item Average Equal Weighted: 1932-07, 1932-08, 1933-04, 1933-05, and 1939-09
  \item Average Value Weighted: 1932-07, 1932-08, 1933-04, and 1933-05
\end{itemize} 
\\

The influential observations are noticeable from the plots of momentum factor MOM against portfolio returns in appendix F. The influential observations have the largest difference between the portfolio return and the MOM factor, in which the corresponding portfolio return are always the peaks and the corresponding MOM factor are always the troughs. For instance, as shown in figure \ref{fig:f2.1} for portfolio 2 with average equal weighted method, the five influential observations are noticeable by inspecting the graphical distance between the MOM factor and the portfolio return.
\\

It was not a coincidence that the influential observations happened in those five months. In fact, they were linked to the major historical events \cite{hist1}. 
\begin{itemize}
  \item \textbf{July - August 1932:} On July 8th, 1932, the Dow reached its lowest level of the 20th century and did not return to pre-1929 levels until the 23rd November, 1954.
  \item \textbf{April - May 1933:} FDR took office in March 1933. In the first 100 days many acts and agencies were introduced  which were to form the basis of the New Deal. The acts and agencies included but not limited to the Economy Act, Emergency Banking Act, National Industrial Recovery Act, Agricultural Adjustment Act, and FDIC. The New Deal made huge impact on economy in years to come.
  \item \textbf{September 1939:} The start of World War II. 
\end{itemize} 
\\

In addition, momentum was noticeable on Dow Jones before and after July 1932. There was tendency for falling asset prices to keep falling from 1929 to July 1932 and tendency for rising asset prices to keep rising from July 1932 to 1937. 
\\

Furthermore, momentum was noticeable on FT 30 but not Dow Jones before and after September 1939. There was tendency for falling asset prices to keep falling from 1936 to September 1939 and tendency for rising asset prices to keep rising from September 1939 to 1946. 
\\

In the next chapter, we will compare the the $\beta$ estimates and the t-statistics between the full dataset and the same dataset with the set of influential observations excluded. 

\begin{figure}[h!] 
\caption{Plot of MOM vs Monthly Portfolio Return (e.g. AEW Portfolio 2)}
\label{fig:f2.3}
<<echo=FALSE>>=
momreturnplotFun(data.l, 1, 2)
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Influence of Influential Observation}

In this chapter, we use the influence measures of Cook's distance and covariance ratio to measure the influence of an observation on the $\beta$ estimates and the t-statistics respectively. We also use both influence measures to examine the impact on the OLS assumptions. To do so, we compare the results of the multiple linear regression on the full dataset and the same dataset with the set of influential observations excluded. \textbf{Panel A} refers to the regression result on the full dataset, \textbf{Panel B} refers to the regression result on the full dataset with the influential points identified by Cook's distance excluded, and \textbf{Panel C} refers to the regression result on the full dataset with the influential points identified by covariance ratio excluded.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Influence on Beta Estimates}

Table \ref{table:t3.1} summarizes the key findings on each of the $\beta$ estimates for both weighting methods. Figures \ref{fig:f3.1} and \ref{fig:f3.2} show the plots of each $\beta$ estimate against the portfolio rankings for the average equal weighted method and the average value weighted method respectively. In sum, the $\beta_{RMRF}$ and $\beta_{SMB}$ estimates are positive and convex on the porfolio rankings for both panels, whereas the $\beta_{HML}$ estimates are positive and concave on the porfolio rankings for both panels. As we should have expected, the $\beta_{MOM}$ estimates are monotonically increasing on the portfolio rankings because the portfolios are constructed monthly using NYSE prior (2-12) return decile breakpoints in which stocks with the lowest NYSE prior (2-12) return comprise decile 1 (or portfolio 1) and stocks with the highest comprise decile 10 (or portfolio 10). Note that the $\beta_{MOM}$ estimates are positive for the average equal weighted method, but for the average value weighted method, the $\beta_{MOM}$ estimates are positive for portfolio ranking less than 7 and negative otherwise.
\\

By excluding the influential observations identified by Cook's distance, the $\beta_{RMRF}$ and $\beta_{SMB}$ estimates are more convex on the porfolio rankings for the average equal weighted method. In addition, $\beta_{HML}$ estimates are more concave on the porfolio rankings for the average equal weighted method. For the rest of the $\beta$ estimates, the graphical relationships and the percentage changes are insignificant. Table \ref{table:t3.2} summarizes the comparisons and the bottom right plots in figures \ref{fig:f3.1} and \ref{fig:f3.2} show the influence in terms of percentage changes.

\begin{table} 
\caption{$\hat{\beta}$ vs Portfolio Ranking}
\label{table:t3.1}
\begin{center}
\begin{tabular}{|c|p{6cm}|p{6cm}|}
\hline
\textbf{$\hat{\beta}$} & \textbf{AEW} & \textbf{AVW} \\
\hline
{$\beta_{RMRF}$} & Positive; Convex on portfolio ranking & Positive; Convex on portfolio ranking \\
\hline
{$\beta_{SMB}$} & Positive; Convex on portfolio ranking & Positive; Convex on portfolio ranking \\
\hline
{$\beta_{HML}$} & Positive; Concave on portfolio ranking & Positive; Concave on portfolio ranking \\
\hline
{$\beta_{MOM}$} & Positive; Monotonically increasing on portfolio ranking & Positive for portfolio ranking less than 7 and negative otherwise; Monotonically increasing on portfolio ranking \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{table} 
\caption{$\hat{\beta}$ vs Portfolio Ranking: Change from Panel A to Panel B}
\label{table:t3.2}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\textbf{$\hat{\beta}$} & \textbf{AEW} & \textbf{AVW} \\
\hline
$\beta_{RMRF}$ & More convex & Insignificant change \\
\hline
$\beta_{SMB}$ & More convex & Insignificant change \\
\hline
$\beta_{HML}$ & More concave & Insignificant change \\
\hline
$\beta_{MOM}$ & Insignificant change & Insignificant change \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}[h!]
\caption{(AEW) Plots of $\beta$}
\label{fig:f3.1}
<<echo=FALSE>>=
par(mfrow=c(3,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta.Ae[,2], beta.Be[,2])
max_y <- max(beta.Ae[,2], beta.Be[,2])
plot(beta.Ae[,1], beta.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Intercept",
     main = "(AEW) Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Ae[,3], beta.Be[,3])
max_y <- max(beta.Ae[,3], beta.Be[,3])
plot(beta.Ae[,1], beta.Ae[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(RMRF)",
     main = "(AEW) Beta of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Ae[,4], beta.Be[,4])
max_y <- max(beta.Ae[,4], beta.Be[,4])
plot(beta.Ae[,1], beta.Ae[,4], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(SMB)",
     main = "(AEW) Beta of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Ae[,5], beta.Be[,5])
max_y <- max(beta.Ae[,5], beta.Be[,5])
plot(beta.Ae[,1], beta.Ae[,5], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(HML)",
     main = "(AEW) Beta of HML")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Ae[,6], beta.Be[,6])
max_y <- max(beta.Ae[,6], beta.Be[,6])
plot(beta.Ae[,1], beta.Ae[,6], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(MOM)",
     main = "(AEW) Beta of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.ABe[,2:6])
max_y <- max(beta.ABe[,2:6])
plot(beta.ABe[,1], beta.ABe[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AEW) % Change from Panel A to Panel B")
axis(1, at=1:10, lab=1:10)
lines(beta.ABe[,1], beta.ABe[,3], type="o", col="green")
lines(beta.ABe[,1], beta.ABe[,4], type="o", col="blue")
lines(beta.ABe[,1], beta.ABe[,5], type="o", col="purple")
lines(beta.ABe[,1], beta.ABe[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y, 
       c("Intercept", "Beta(RMRF)", "Beta(SMB)", "Beta(HML)", "Beta(MOM)"), 
       cex=0.6, col=c("red", "green", "blue", "purple", "brown"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\caption{(AVW) Plots of $\beta$}
\label{fig:f3.2}
<<echo=FALSE>>=
par(mfrow=c(3,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta.Av[,2], beta.Bv[,2])
max_y <- max(beta.Av[,2], beta.Bv[,2])
plot(beta.Av[,1], beta.Av[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Intercept",
     main = "(AVW) Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Av[,3], beta.Bv[,3])
max_y <- max(beta.Av[,3], beta.Bv[,3])
plot(beta.Av[,1], beta.Av[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(RMRF)",
     main = "(AVW) Beta of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Av[,4], beta.Bv[,4])
max_y <- max(beta.Av[,4], beta.Bv[,4])
plot(beta.Av[,1], beta.Av[,4], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(SMB)",
     main = "(AVW) Beta of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Av[,5], beta.Bv[,5])
max_y <- max(beta.Av[,5], beta.Bv[,5])
plot(beta.Av[,1], beta.Av[,5], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(HML)",
     main = "(AVW) Beta of HML")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.Av[,6], beta.Bv[,6])
max_y <- max(beta.Av[,6], beta.Bv[,6])
plot(beta.Av[,1], beta.Av[,6], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(MOM)",
     main = "(AVW) Beta of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta.ABv[,2:6])
max_y <- max(beta.ABv[,2:6])
plot(beta.ABv[,1], beta.ABv[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AVW) % Change from Panel A to Panel B")
axis(1, at=1:10, lab=1:10)
lines(beta.ABv[,1], beta.ABv[,3], type="o", col="green")
lines(beta.ABv[,1], beta.ABv[,4], type="o", col="blue")
lines(beta.ABv[,1], beta.ABv[,5], type="o", col="purple")
lines(beta.ABv[,1], beta.ABv[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y, 
       c("Intercept", "Beta(RMRF)", "Beta(SMB)", "Beta(HML)", "Beta(MOM)"), 
       cex=0.6, col=c("red", "green", "blue", "purple", "brown"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\section{Influence on t-statistics}

In this section, we want to examine the influence on the t-statistic and also the Newey-West t-statistic. As we have discussed in section 1.6, we use the Newey-West t-statistic as defined in equation \eqref{eq:eq1.17} to overcome autocorrelation and heteroskedasticity in the random errors in the Fama French Carhart Four Model \eqref{eq:eq1.1}. Tables \ref{table:t3.3} and \ref{table:t3.4} summarize the key findings on both versions of t-statistic of the $\beta$ estimates for the average equal weighted method and the average value weighted method respectively. Figures \ref{fig:f3.4}, \ref{fig:f3.5}, \ref{fig:f3.6}, and \ref{fig:f3.7} show the plots of both versions of t-statistic against the portfolio rankings for $\beta_{RMRF}$, $\beta_{SMB}$, $\beta_{HML}$, and $\beta_{MOM}$ respectively. I also include the same kind of plots for the intercept in figure \ref{fig:f3.3}.
\\

For the average equal weighted method, both versions of t-statistic are positive and highly significant for all the $\beta$ estimates except for the $\beta_{MOM}$ estimates in portfolios 7-10. As for the average value weighted method, both versions of t-statistic are positive for all the $\beta$ estimates except for the $\beta_{SMB}$ estimates in portfolios 3-9 and the $\beta_{MOM}$ estimates in portfolios 7-10. The $\beta_{SMB}$ estimates are highly significant except for portfolios 3 and 9, whereas the $\beta_{HML}$ estimates are highly significant except for portfolios 1 and 10. In addition, $\beta_{MOM}$ estimates are monotonically increasing on portfolio ranking starting from portfolios 2 or 3 across the board.
\\

By excluding the influential observations identified by covariance ratio, the sign of both t-statistic, the level of statistical significance, and the graphical relationship between the t-statistic and the porfolio rankings remain the same at large. However, when we look into the percentage change as shown in figure \ref{fig:f3.8}, the percentage changes from panel A to panel C are much higher for the Newey-West t-statistic than the ones for t-statistic across the board.

\begin{table} 
\caption{(AEW) t-statistic and Newey-West t-statistic for $\hat{\beta}$}
\label{table:t3.3}
\begin{center}
\begin{tabular}{|c|p{7cm}|p{7cm}|}
\hline
\textbf{$\hat{\beta}$} & \textbf{t-statistic} & \textbf{Newey-West t-statistic} \\
\hline
{$\beta_{RMRF}$} & Positive; Highly significant; Concave on portfolio ranking & Positive; Highly significant; Multimodal on portfolio ranking \\
\hline
{$\beta_{SMB}$} & Positive; Highly significant; Multimodal on portfolio ranking & Positive; Highly significant; Multimodal on portfolio ranking \\
\hline
{$\beta_{HML}$} & Positive; Highly significant; Near concave on portfolio ranking & Positive; Highly significant; Near concave on portfolio ranking \\
\hline
{$\beta_{MOM}$} & Positive for portfolio ranking less than 7 and negative otherwise; Highly significant; Monotonically increasing on portfolio ranking starting from portfolios 2 or 3 &  Positive for portfolio ranking less than 7 and negative otherwise; Highly significant; Monotonically increasing on portfolio ranking starting from portfolios 2 or 3 \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{table} 
\caption{(AVW) t-statistic and Newey-West t-statistic for $\hat{\beta}$}
\label{table:t3.4}
\begin{center}
\begin{tabular}{|c|p{7cm}|p{7cm}|}
\hline
\textbf{$\hat{\beta}$} & \textbf{t-statistic} & \textbf{Newey-West t-statistic} \\
\hline
{$\beta_{RMRF}$} & Positive; Highly significant; Multimodal on portfolio ranking & Positive; Highly significant; Multimodal on portfolio ranking \\
\hline
{$\beta_{SMB}$} & Positive for portfolios 1, 2, 10, and negative otherwise; Highly significant except for portfolios 3 and 9; Convex on portfolio ranking & Positive for portfolios 1, 2, 10, and negative otherwise; Highly significant except for portfolios 3 and 9; Convex on portfolio ranking \\
\hline
{$\beta_{HML}$} & Positive; Highly significant except for portfolios 1 and 10; Near concave on portfolio ranking & Positive; Highly significant except for portfolios 1 and 10; Near concave on portfolio ranking \\
\hline
{$\beta_{MOM}$} & Positive for portfolio ranking less than 7 and negative otherwise; Highly significant; Monotonically increasing on portfolio ranking starting from portfolios 2 or 3 &  Positive for portfolio ranking less than 7 and negative otherwise; Highly significant; Monotonically increasing on portfolio ranking starting from portfolios 2 or 3 \\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}[h!]
\caption{Plots of Intercept t-statistic}
\label{fig:f3.3}
<<echo=FALSE>>=
par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,2], beta_t.Ce[,2])
max_y <- max(beta_t.Ae[,2], beta_t.Ce[,2])
plot(beta_t.Ae[,1], beta_t.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(Intercept)",
     main = "(AEW) T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_t.Av[,2], beta_t.Cv[,2])
max_y <- max(beta_t.Av[,2], beta_t.Cv[,2])
plot(beta_t.Av[,1], beta_t.Av[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(Intercept)",
     main = "(AVW) T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,2], beta_nwt.Ce[,2])
max_y <- max(beta_nwt.Ae[,2], beta_nwt.Ce[,2])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(Intercept)",
     main = "(AEW) NW T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Av[,2], beta_nwt.Cv[,2])
max_y <- max(beta_nwt.Av[,2], beta_nwt.Cv[,2])
plot(beta_nwt.Av[,1], beta_nwt.Av[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(Intercept)",
     main = "(AVW) NW T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\caption{Plots of $\beta_{RMRF}$ t-statistic}
\label{fig:f3.4}
<<echo=FALSE>>=
par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,3], beta_t.Ce[,3])
max_y <- max(beta_t.Ae[,3], beta_t.Ce[,3])
plot(beta_t.Ae[,1], beta_t.Ae[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(RMRF)",
     main = "(AEW) T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_t.Av[,3], beta_t.Cv[,3])
max_y <- max(beta_t.Av[,3], beta_t.Cv[,3])
plot(beta_t.Av[,1], beta_t.Av[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(RMRF)",
     main = "(AVW) T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,3], beta_nwt.Ce[,3])
max_y <- max(beta_nwt.Ae[,3], beta_nwt.Ce[,3])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(RMRF)",
     main = "(AEW) NW T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Av[,3], beta_nwt.Cv[,3])
max_y <- max(beta_nwt.Av[,3], beta_nwt.Cv[,3])
plot(beta_nwt.Av[,1], beta_nwt.Av[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(RMRF)",
     main = "(AVW) NW T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\caption{Plots of $\beta_{SMB}$ t-statistic}
\label{fig:f3.5}
<<echo=FALSE>>=
par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,4], beta_t.Ce[,4])
max_y <- max(beta_t.Ae[,4], beta_t.Ce[,4])
plot(beta_t.Ae[,1], beta_t.Ae[,4], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(SMB)",
     main = "(AEW) T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_t.Av[,4], beta_t.Cv[,4])
max_y <- max(beta_t.Av[,4], beta_t.Cv[,4])
plot(beta_t.Av[,1], beta_t.Av[,4], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(SMB)",
     main = "(AVW) T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,4], beta_nwt.Ce[,4])
max_y <- max(beta_nwt.Ae[,4], beta_nwt.Ce[,4])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,4], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(SMB)",
     main = "(AEW) NW T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Av[,4], beta_nwt.Cv[,4])
max_y <- max(beta_nwt.Av[,4], beta_nwt.Cv[,4])
plot(beta_nwt.Av[,1], beta_nwt.Av[,4], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(SMB)",
     main = "(AVW) NW T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\caption{Plots of $\beta_{HML}$ t-statistic}
\label{fig:f3.6}
<<echo=FALSE>>=
par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,5], beta_t.Ce[,5])
max_y <- max(beta_t.Ae[,5], beta_t.Ce[,5])
plot(beta_t.Ae[,1], beta_t.Ae[,5], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(HML)",
     main = "(AEW) T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_t.Av[,5], beta_t.Cv[,5])
max_y <- max(beta_t.Av[,5], beta_t.Cv[,5])
plot(beta_t.Av[,1], beta_t.Av[,5], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(HML)",
     main = "(AVW) T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,5], beta_nwt.Ce[,5])
max_y <- max(beta_nwt.Ae[,5], beta_nwt.Ce[,5])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,5], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(HML)",
     main = "(AEW) NW T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Av[,5], beta_nwt.Cv[,5])
max_y <- max(beta_nwt.Av[,5], beta_nwt.Cv[,5])
plot(beta_nwt.Av[,1], beta_nwt.Av[,5], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(HML)",
     main = "(AVW) NW T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\caption{Plots of $\beta_{MOM}$ t-statistic}
\label{fig:f3.7}
<<echo=FALSE>>=
par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,6], beta_t.Ce[,6])
max_y <- max(beta_t.Ae[,6], beta_t.Ce[,6])
plot(beta_t.Ae[,1], beta_t.Ae[,6], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(MOM)",
     main = "(AEW) T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_t.Av[,6], beta_t.Cv[,6])
max_y <- max(beta_t.Av[,6], beta_t.Cv[,6])
plot(beta_t.Av[,1], beta_t.Av[,6], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(MOM)",
     main = "(AVW) T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,6], beta_nwt.Ce[,6])
max_y <- max(beta_nwt.Ae[,6], beta_nwt.Ce[,6])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,6], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(MOM)",
     main = "(AEW) NW T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.Av[,6], beta_nwt.Cv[,6])
max_y <- max(beta_nwt.Av[,6], beta_nwt.Cv[,6])
plot(beta_nwt.Av[,1], beta_nwt.Av[,6], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(MOM)",
     main = "(AVW) NW T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[h!]
\caption{Plots of Percentage Change in t-statistic}
\label{fig:f3.8}
<<echo=FALSE>>=
par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.ACe[,2:6])
max_y <- max(beta_t.ACe[,2:6])
plot(beta_t.ACe[,1], beta_t.ACe[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AEW) % Change in T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_t.ACe[,1], beta_t.ACe[,3], type="o", col="green")
lines(beta_t.ACe[,1], beta_t.ACe[,4], type="o", col="blue")
lines(beta_t.ACe[,1], beta_t.ACe[,5], type="o", col="purple")
lines(beta_t.ACe[,1], beta_t.ACe[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y, 
       c("Intercept", "T(RMRF)", "T(SMB)", "T(HML)", "T(MOM)"), 
       cex=0.6, col=c("red", "green", "blue", "purple", "brown"), 
       lty=1, bty="n")

min_y <- min(beta_t.ACv[,2:6])
max_y <- max(beta_t.ACv[,2:6])
plot(beta_t.ACv[,1], beta_t.ACv[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AVW) % Change in T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_t.ACv[,1], beta_t.ACv[,3], type="o", col="green")
lines(beta_t.ACv[,1], beta_t.ACv[,4], type="o", col="blue")
lines(beta_t.ACv[,1], beta_t.ACv[,5], type="o", col="purple")
lines(beta_t.ACv[,1], beta_t.ACv[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y, 
       c("Intercept", "T(RMRF)", "T(SMB)", "T(HML)", "T(MOM)"), 
       cex=0.6, col=c("red", "green", "blue", "purple", "brown"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.ACe[,2:6])
max_y <- max(beta_nwt.ACe[,2:6])
plot(beta_nwt.ACe[,1], beta_nwt.ACe[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AEW) % Change in NW T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,3], type="o", col="green")
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,4], type="o", col="blue")
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,5], type="o", col="purple")
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y, 
       c("Intercept", "NWT(RMRF)", "NWT(SMB)", "NWT(HML)", "NWT(MOM)"), 
       cex=0.6, col=c("red", "green", "blue", "purple", "brown"), 
       lty=1, bty="n")

min_y <- min(beta_nwt.ACv[,2:6])
max_y <- max(beta_nwt.ACv[,2:6])
plot(beta_nwt.ACv[,1], beta_nwt.ACv[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AVW) % Change in NW T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,3], type="o", col="green")
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,4], type="o", col="blue")
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,5], type="o", col="purple")
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y, 
       c("Intercept", "NWT(RMRF)", "NWT(SMB)", "NWT(HML)", "NWT(MOM)"), 
       cex=0.6, col=c("red", "green", "blue", "purple", "brown"), 
       lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\section{Influence on the OLS Assumptions}
In this section, we want to examine the impact of the influential observations on the OLS assumptions, in particular, normality and variance of residuals.

\subsection{Lilliefors Test for Normality}
We apply the Lilliefors test on panels B and C to test the null hypothesis that the residuals come from the Gaussian distribution. Tables \ref{table:t3.5} and \ref{table:t3.6} show the test results for panels B and C respectively. The Lilliefors test results for panel B remain the same as the ones for panel A. However, for panel C, the normality assumption seem to be valid for 3 and 8 out of 10 portfolios for the average equal weighted and average value weighted methods respectively.
\\

<<echo=FALSE, results="asis">>=
print(
  xtable(normtestpval.B, 
         caption = "(Panel B) Lilliefors Test for Normality with 0.05 Alpha Level",
         label = "table:t3.5", display = c(rep('d', 2), rep('e', 2), rep('s', 2))),
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

<<echo=FALSE, results="asis">>=
print(
  xtable(normtestpval.C, 
         caption = "(Panel C) Lilliefors Test for Normality with 0.05 Alpha Level",
         label = "table:t3.6", display = c(rep('d', 2), rep('e', 2), rep('s', 2))),
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

\subsection{Non-constant Variance Score Test}
We apply the non-constant variance score test on panels B and C to test the null hypothesis that the residuals are distributed with constant variance. Tables \ref{table:t3.7} and \ref{table:t3.8} show the test results for panels B and C respectively. For panel B, there is 1 additional portfolio seem to be distributed with constant variance for each weighting methods. Whereas for panel C, there are 7 and 5 additional portfolios seem to be distributed with constant variance for the average equal weighted and average value weighted methods respectively.
\\

<<echo=FALSE, results="asis">>=
print(
  xtable(ncvtestpval.B, 
         caption = "(Panel B) Non-constant Variance Score Test with 0.05 Alpha Level",
         label = "table:t3.7", display = c(rep('d', 2), rep('e', 2), rep('s', 2))),
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

<<echo=FALSE, results="asis">>=
print(
  xtable(ncvtestpval.C, 
         caption = "(Panel C) Non-constant Variance Score Test with 0.05 Alpha Level",
         label = "table:t3.8", display = c(rep('d', 2), rep('e', 2), rep('s', 2))),
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Portfolio Optimization}

As a continuation from the last chapter, we want to examine the impact of the influence observations on the expected monthly return and the standard deviation of the monthly return by comparing their results among panels A, B, and C. Then, with the two different weighting methods for consideration, our logical next step is to compare their investment performances by looking at their mean-variance optimal portfolios.  

\section{Mean and SD of Monthly Portfolio Returns}

Figure \ref{fig:f4.1} shows the plots of the expected monthly portfolio returns, their standard deviations, and the percentage change from panel A to panels B and C against the portfolio rankings for both weighting methods. In sum, for both weighting methods,
\begin{itemize}
  \item The expected monthly portfolio return increases as the portfolio ranking increases. 
  \item The standard deviations of monthly portfolio returns are convex on the portfolio ranking.  
  \item The percentage changes for the expected monthly portfolio returns and their standard deviations are negative from panel A to panels B and C, except for the expected monthly portfolio returns for portfolios 9 and 10 with the average value weighted method in panel C.
  \item The absolute percentage changes from panel A to panels B and C for the two estimates (mean and SD) are diminishing on the portfolio rankings at large. 
\end{itemize}
\\

To go one step further, we want to look into the relationship between the expected monthly portfolio returns and their standard deviations for each of the weighting methods and panels. As shown in the top two plots in figure \ref{fig:f4.2}, for both weighting methods, the risk and return curve shifts to the left when we change from panel A to panel B. Similarly, for both weighting methods, the risk and return curve shifts to the left when we change from either panel A or panel B to panel C. These results are what we have expected because the influential observations have abnormal high monthly portfolio returns and therefore their removal would reduce the risk and return. The bottom two plots show that the risk and return curve shifts to the left when we change from the average equal weighted method to the average value weighted method in all panels. The explanation will be given in the latter section \textit{Optimal Portfolio Weights}.
\\

\begin{figure}[h!]
\caption{Plots of Mean and SD of Monthly Portfolio Returns}
\label{fig:f4.1}
<<echo=FALSE>>=
### Fill by column
par(mfcol=c(3,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,3))

min_y <- min(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
max_y <- max(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
plot(meansd.Ae[,1], meansd.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Mean",
     main = "(AEW) Mean of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Be[,1], meansd.Be[,2], type="o", col="red")
lines(meansd.Ce[,1], meansd.Ce[,2], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6, 
       col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
max_y <- max(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
plot(meansd.Ae[,1], meansd.Ae[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "SD",
     main = "(AEW) SD of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Be[,1], meansd.Be[,3], type="o", col="red")
lines(meansd.Ce[,1], meansd.Ce[,3], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6, 
       col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.ABe[,2:3], meansd.ACe[,2:3])
max_y <- max(meansd.ABe[,2:3], meansd.ACe[,2:3])
plot(meansd.ABe[,1], meansd.ABe[,2], type = "o", col = "skyblue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AEW) % Change from Panel A to Panels B & C")
axis(1, at=1:10, lab=1:10)
lines(meansd.ACe[,1], meansd.ACe[,2], type="o", col="orange")
lines(meansd.ABe[,1], meansd.ABe[,3], type="o", col="purple")
lines(meansd.ACe[,1], meansd.ACe[,3], type="o", col="brown")
legend(10.5, max_y, c("Mean (B)", "Mean (C)", "SD (B)", "SD (C)"), cex=0.6, 
       col=c("skyblue", "orange", "purple", "brown"), lty=1, bty="n")

min_y <- min(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
max_y <- max(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
plot(meansd.Av[,1], meansd.Av[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Mean",
     main = "(AVW) Mean of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Bv[,1], meansd.Bv[,2], type="o", col="red")
lines(meansd.Cv[,1], meansd.Cv[,2], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6, 
       col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
max_y <- max(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
plot(meansd.Av[,1], meansd.Av[,3], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "SD",
     main = "(AVW) SD of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Bv[,1], meansd.Bv[,3], type="o", col="red")
lines(meansd.Cv[,1], meansd.Cv[,3], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6, 
       col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.ABv[,2:3], meansd.ACv[,2:3])
max_y <- max(meansd.ABv[,2:3], meansd.ACv[,2:3])
plot(meansd.ABv[,1], meansd.ABv[,2], type = "o", col = "skyblue", 
     ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
     main = "(AVW) % Change from Panel A to Panels B & C")
axis(1, at=1:10, lab=1:10)
lines(meansd.ACv[,1], meansd.ACv[,2], type="o", col="orange")
lines(meansd.ABv[,1], meansd.ABv[,3], type="o", col="purple")
lines(meansd.ACv[,1], meansd.ACv[,3], type="o", col="brown")
legend(10.5, max_y, c("Mean (B)", "Mean (C)", "SD (B)", "SD (C)"), cex=0.6, 
       col=c("skyblue", "orange", "purple", "brown"), lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))
@
\end{figure}

\begin{figure}[h!]
\caption{Plots of Mean versus SD of Monthly Portfolio Returns}
\label{fig:f4.2}
<<echo=FALSE>>=
par(mfrow=c(3,2))

min_x <- min(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
max_x <- max(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
min_y <- min(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
max_y <- max(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
plot(meansd.Ae[,3], meansd.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlim = c(min_x,max_x),
     xlab = "SD", ylab = "Mean", main = "(AEW) Mean vs SD")
lines(meansd.Be[,3], meansd.Be[,2], type="o", col="red")
lines(meansd.Ce[,3], meansd.Ce[,2], type="o", col="brown")
legend("topright", c("Panel A", "Panel B", "Panel C"), cex=0.6, 
       col=c("blue", "red", "brown"), lty=1, bty="n")

min_x <- min(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
max_x <- max(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
min_y <- min(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
max_y <- max(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
plot(meansd.Av[,3], meansd.Av[,2], type = "o", col = "green", 
     ylim = c(min_y,max_y), xlim = c(min_x,max_x),
     xlab = "SD", ylab = "Mean", main = "(AVW) Mean vs SD")
lines(meansd.Bv[,3], meansd.Bv[,2], type="o", col="purple")
lines(meansd.Cv[,3], meansd.Cv[,2], type="o", col="orange")
legend("topright", c("Panel A", "Panel B", "Panel C"), cex=0.6, 
       col=c("green", "purple", "orange"), lty=1, bty="n")

min_x <- min(meansd.Ae[,3], meansd.Av[,3])
max_x <- max(meansd.Ae[,3], meansd.Av[,3])
min_y <- min(meansd.Ae[,2], meansd.Av[,2])
max_y <- max(meansd.Ae[,2], meansd.Av[,2])
plot(meansd.Ae[,3], meansd.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlim = c(min_x,max_x),
     xlab = "SD", ylab = "Mean", main = "(Panel A) Mean vs SD")
lines(meansd.Av[,3], meansd.Av[,2], type="o", col="green")
legend("topright", c("AEW", "AVW"), cex=0.6, col=c("blue", "green"), 
       lty=1, bty="n")

min_x <- min(meansd.Be[,3], meansd.Bv[,3])
max_x <- max(meansd.Be[,3], meansd.Bv[,3])
min_y <- min(meansd.Be[,2], meansd.Bv[,2])
max_y <- max(meansd.Be[,2], meansd.Bv[,2])
plot(meansd.Be[,3], meansd.Be[,2], type = "o", col = "red", 
     ylim = c(min_y,max_y), xlim = c(min_x,max_x),
     xlab = "SD", ylab = "Mean", main = "(Panel B) Mean vs SD")
lines(meansd.Bv[,3], meansd.Bv[,2], type="o", col="purple")
legend("topright", c("AEW", "AVW"), cex=0.6, col=c("red", "purple"), 
       lty=1, bty="n")

min_x <- min(meansd.Ce[,3], meansd.Cv[,3])
max_x <- max(meansd.Ce[,3], meansd.Cv[,3])
min_y <- min(meansd.Ce[,2], meansd.Cv[,2])
max_y <- max(meansd.Ce[,2], meansd.Cv[,2])
plot(meansd.Ce[,3], meansd.Ce[,2], type = "o", col = "brown", 
     ylim = c(min_y,max_y), xlim = c(min_x,max_x),
     xlab = "SD", ylab = "Mean", main = "(Panel C) Mean vs SD")
lines(meansd.Cv[,3], meansd.Cv[,2], type="o", col="orange")
legend("topright", c("AEW", "AVW"), cex=0.6, col=c("brown", "orange"), 
       lty=1, bty="n")

min_x <- min(meansd.Ae[,3], meansd.Av[,3], 
             meansd.Be[,3], meansd.Bv[,3],
             meansd.Ce[,3], meansd.Cv[,3])
max_x <- max(meansd.Ae[,3], meansd.Av[,3], 
             meansd.Be[,3], meansd.Bv[,3],
             meansd.Ce[,3], meansd.Cv[,3])
min_y <- min(meansd.Ae[,2], meansd.Av[,2], 
             meansd.Be[,2], meansd.Bv[,2],
             meansd.Ce[,2], meansd.Cv[,2])
max_y <- max(meansd.Ae[,2], meansd.Av[,2], 
             meansd.Be[,2], meansd.Bv[,2],
             meansd.Ce[,2], meansd.Cv[,2])
plot(meansd.Ae[,3], meansd.Ae[,2], type = "o", col = "blue", 
     ylim = c(min_y,max_y), xlim = c(min_x,max_x),
     xlab = "SD", ylab = "Mean", main = "Mean vs SD")
lines(meansd.Av[,3], meansd.Av[,2], type="o", col="green")
lines(meansd.Be[,3], meansd.Be[,2], type = "o", col = "red")
lines(meansd.Bv[,3], meansd.Bv[,2], type="o", col="purple")
lines(meansd.Ce[,3], meansd.Ce[,2], type = "o", col = "brown")
lines(meansd.Cv[,3], meansd.Cv[,2], type="o", col="orange")
legend("topright", 
       c("AEW (A)", "AVW (A)", "AEW (B)", "AVW (B)", "AEW (C)", "AVW (C)"), 
       cex=0.6, 
       col=c("blue", "green", "red", "purple", "brown", "orange"), 
       lty=1, bty="n")

par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\section{Mean-Variance Optimal Portfolio}

In this section, we want to compare the investment performances of two mean-variance optimal portfolios, in which one portfolio is comprised of the ten momentum portfolios formed with the average equal weighted method and the other portfolio is comprised of the ten momentum portfolios formed with the average value weighted method. 
\\

To find a mean-variance optimal portfolio with ten assets \cite{portfolio1}, i.e. the ten momentum portfolios, first, we compute the mean rates of return $\bar{r}_1, \bar{r}_2, \ldots, \bar{r}_{10}$ and the covariances $\sigma_{ij}$, $i, j = 1, 2, \ldots, 10$ for the ten momentum portfolios, using all the months in our dataset. Our portfolio is defined by a set of ten weights $w_{i}$, $i = 1, 2, \ldots, 10$, that sum to 1 (Negative weight corresponds to short selling). Then, we fix the mean value at some arbitrary value $\bar{r}$ and find the feasible portfolio of minimum variance that has this mean. The mathematical problem of solving for the mean-variance optimal portfolio as known as the Markowitz problem is equivalent to the following problem:

\begin{center}
$minimize \, \frac{1}{2}\sum_{i,j = 1}^{10}w_{i}w_{j}\sigma_{ij}$
\end{center}
\\

\begin{center}
$subject \, to \, \sum_{i = 1}^{10}w_{i}{\bar{r}}_{i} = \bar{r}$
\end{center}
\\

\begin{center}
$\sum_{i = 1}^{10}w_{i} = 1$
\end{center}
\\

To find the solutions for the Markowitz problem, we form the \textit{Lagrangian} $L$ with the \textit{Lagrange multipliers} $\lambda$ and $\mu$. 
\begin{equation} 
\label{eq:eq4.1}
L = \frac{1}{2}\sum_{i,j = 1}^{10}w_{i}w_{j}\sigma_{ij} - \lambda\left(\sum_{i = 1}^{10}w_{i}{\bar{r}}_{i} - \bar{r}\right) - \mu\left(\sum_{i = 1}^{10}w_{i} - 1\right)
\end{equation}
\\

After setting the first derivatives of the \textit{Lagrangian} $L$ with respect to $\lambda$, $\mu$, and the 10 portfolio weights $w_i$ ($i = 1, 2, \ldots, 10$) equal to zero, we obtain
\begin{equation} 
\label{eq:eq4.2}
\sum_{j = 1}^{10}\sigma_{ij}w_{j} - \lambda{\bar{r}}_{i} - \mu = 0 \; for \, i = 1, 2, \ldots, 10
\end{equation}

\begin{equation} 
\label{eq:eq4.3}
\sum_{i = 1}^{10}w_{i}{\bar{r}}_{i} = \bar{r}
\end{equation}

\begin{equation} 
\label{eq:eq4.4}
\sum_{i = 1}^{10}w_{i} = 1
\end{equation}
\\

With equations \eqref{eq:eq4.2}, \eqref{eq:eq4.3}, and \eqref{eq:eq4.4}, we can generate the mean-variance optimal set of portfolios as shown in figure \ref{fig:f4.3}. The top two plots compares the mean-variance optimal sets with selling short allowed and going long only, in which the plot on the top refers to the average equal weighted method and the plot in the middle refers to the average value weighted method. Both plots show that the mean-variance optimal set shifts to the left when selling short is allowed, which is what we have expected because the feasible solution set for going long only is the subset of the feasible solution set for selling short allowed.
\\

Absent risk-free borrowing, the efficient frontier is the upper arc of the parabola (the mean-variance optimal set) because for any portfolio on the lower arc you can always find a portfolio on the upper arc with the same standard deviation and higher expected return. Once a risk-free rate is introduced, the efficient frontier is a line that is tangent to the upper arc of the parabola. \cite{portfolio1}.
\\

With the assumption of selling short allowed, the plot at the bottom compares the mean-variance optimal sets for both weighting method. The plot shows that the mean-variance optimal set shifts to the left when we change from the average equal weighted method to the average value weighted method, but with the intersection near the expected monthly return of $0.018$ and its standard deviation of $0.060$. Therefore, if our risk appetite (or standard deviation) is lower than $0.060$, we should hold the optimal portfolio formed with the average value weighted method, otherwise, we should hold the optimal portfolio formed with the average equal weighted method.
\\

The result that the mean-variance optimal set shifts to the left when we change from the average equal weighted method to the average value weighted method can be explained by Banz's finding \cite{weight1} that smaller firms have had higher risk adjusted returns, on average, than larger firms. The average equal weighted method put much more weight on the smaller firms than the larger firms because we have numerous of smaller firms in the market. Therefore, the lowest level of risk that is achievable with the average equal weighted method is higher than the one that is achievable with the average value weighted method. However, with the effect of diversification and selling short allowed, we can achieve higher expected monthly return with risk or standard deviation greater than $0.060$ using the average equal weighted method.
\\

\begin{figure}[h!]
\caption{Plots of Mean-Variance Optimal Portfolio}
\label{fig:f4.3}
<<echo=FALSE>>=
par(mfrow=c(3,1))

### (AEW) Mean-Variance Optimal Portfolio
plot(effpfT.e[,2], effpfT.e[,1], type = "o", col = "red", 
     xlim = c(min_x.e,max_x.e), xlab = "Standard Deviation", 
     ylab = "Expected Monthly Return", 
     main = "(AEW) Mean-Variance Optimal Portfolio")
lines(effpfF.e[,2], effpfF.e[,1], type="o", col="blue")
legend("right", c("Short Allowed", "Long Only"), cex=0.6, 
       col=c("red", "blue"), lty=1, bty="n")

### (AVW) Mean-Variance Optimal Portfolio
plot(effpfT.v[,2], effpfT.v[,1], type = "o", col = "red", 
     xlim = c(min_x.v,max_x.v), xlab = "Standard Deviation", 
     ylab = "Expected Monthly Return", 
     main = "(AVW) Mean-Variance Optimal Portfolio")
lines(effpfF.v[,2], effpfF.v[,1], type="o", col="blue")
legend("right", c("Short Allowed", "Long Only"), cex=0.6, 
       col=c("red", "blue"), lty=1, bty="n")

### (Shorts Allowed) Mean-Variance Optimal Portfolio: AEW vs AVW
plot(effpf.e[,2], effpf.e[,1], type = "o", col = "red", 
     xlim = c(min_x.ev,max_x.ev), xlab = "Standard Deviation", 
     ylab = "Expected Monthly Return", 
     main = "(Shorts Allowed) Mean-Variance Optimal Portfolio: AEW vs AVW")
lines(effpf.v[,2], effpf.v[,1], type="o", col="blue")
legend("right", c("AEW", "AVW"), cex=0.6, 
       col=c("red", "blue"), lty=1, bty="n")

par(mfrow=c(1,1))
@
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

\section{Optimal Portfolio Weights} 

Tables \ref{table:t4.1} and \ref{table:t4.2} contain the solutions for the mean-variance optimal portfolio for the average equal weighted method and the average value weighted method respectively. It is easier to see the relationships among the minimum portfolio set and the corresponding optimal portfolio weights by looking at their correlations. The numerical results for the average equal weighted method and the average value weighted method are shown in tables \ref{table:t4.3} and \ref{table:t4.4} respectively; whereas the graphic results are shown in figures \ref{fig:f4.4} and \ref{fig:f4.5} respectively.
\\

We find that for the average equal weighted method, the optimal weights for portfolios 2, 4, 6, 9, and 10 are positively correlated with the expected month return and risk; whereas for the average value weighted method, the optimal weights for portfolios 2, 4, and 10 are positively correlated with the expected month return and risk. Therefore, the average equal weighted method includes two additional portfolio (9 and 10) for the positive correlation.
\\

<<echo=FALSE, results="asis">>=
print(
  xtable(effpf.e, 
         caption = "(AEW) Mean-Variance Optimal Portfolio: Solution",
         label = "table:t4.1", digits = c(0, rep(3, 2), rep(2, 10)), 
         display = c('e', rep('f', 12))), 
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

<<echo=FALSE, results="asis">>=
print(
  xtable(effpf.v, 
         caption = "(AVW) Mean-Variance Optimal Portfolio: Solution",
         label = "table:t4.2", digits = c(0, rep(3, 2), rep(2, 10)), 
         display = c('e', rep('f', 12))), 
  table.placement = "htbp", caption.placement = "top", 
  include.rownames = FALSE, latex.environment = "center")
@

\begin{figure}[h!]
\caption{(AEW) Scatterplot of the Portfolio's Mean, SD, and Weight}
\label{fig:f4.4}
<<echo=FALSE>>=
pairs(effpfT.e, col="blue")
@
\end{figure}

\begin{figure}[h!]
\caption{(AVW) Scatterplot of the Portfolio's Mean, SD, and Weight}
\label{fig:f4.5}
<<echo=FALSE>>=
pairs(effpfT.v, col="blue")
@
\end{figure}

<<echo=FALSE, results="asis">>=
print(
  xtable(cor(effpf.e), 
         caption = "(AEW) Mean-Variance Optimal Portfolio: Correlation",
         label = "table:t4.3", digits = c(0, rep(1, 12)), 
         display = c('e', rep('f', 12))), 
  table.placement = "htbp", caption.placement = "top", 
  latex.environment = "center")
@

<<echo=FALSE, results="asis">>=
print(
  xtable(cor(effpf.v), 
         caption = "(AVW) Mean-Variance Optimal Portfolio: Correlation",
         label = "table:t4.4", digits = c(0, rep(1, 12)), 
         display = c('e', rep('f', 12))), 
  table.placement = "htbp", caption.placement = "top", 
  latex.environment = "center")
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Conclusion}

Based on the results of Lilliefors test for normality, with $5\%$ $\alpha$ level, the random errors are not Gaussian distributed, therefore the OLS estimates are not equivalent to the maximum-likelihood estimates. For OLS, the random errors do not have to be Gaussian distributed.
\\

Through our regression diagnostic, the random errors appear to be distributed with mean zero. In addition, the random errors appear to be independent of the observed data of the Fama French Carhart four factors. However, based on the times series of residuals, the residuals seem not to be independent and identically distributed because they have larger fluctuations before 1940 and after 1980. The time series plots also suggest that the homoscedasticity (constant variance) assumption is violated, which are confirmed by the results of the non-constant variance score test with $5\%$ alpha level. Our failure of meeting some of the OLS assumptions suggests that our dataset constains influential observations.
\\

Covariance ratio identifies more influential observations than Cook's distance. Out of 1,039 monthly observations, for the average equal weighted method, Cook's distance and covariance ratio identify 5 and 187 influential observations respectively. For the average value weighted method, Cook's distance and covariance ratio identify 4 and 232 influential observations respectively.
\\

We use Cook's distance to examine influence on betas estimates. After excluding the influential observations identified by Cook's distance from the original data set, for the average equal weighted method, the $\beta_{RMRF}$ and $\beta_{SMB}$ estimates are more convex on portfolio rankings and $\beta_{HML}$ estimates are more concave on portfolio rankings. We have also found that all the influential observations identified by Cook's Distance are connected with major historical events. 
\\

We use covariance ratio to examine influence on t-statistics. After excluding the influential observations identified by covariance ratio from the original data set, there are substantial percentage changes for both the Newey-West t-statistics and the standard t-statistics. The Newey-West t-statistics appear to be more sensitive to the elimination of influential observations than the ordinary t-statistics. 
\\

We apply the Lilliefors test on panels B and C to test the null hypothesis that the residuals come from the Gaussian distribution. For panel B, the results of Lilliefors test for normality remain the same as the ones for panel A. However, for panel C, the normality assumption seem to be valid for 3 and 8 out of 10 portfolios for the average equal weighted and average value weighted methods respectively.
\\

We also apply the non-constant variance score test on panels B and C to test the null hypothesis that the residuals are distributed with constant variance. For panel B, there is 1 additional portfolio seem to be distributed with constant variance for each weighting methods. Whereas for panel C, there are 7 and 5 additional portfolios seem to be distributed with constant variance for the average equal weighted and average value weighted methods respectively.
\\

By comparing the two mean-variance optimal portfolios, in which the portfolios are both comprised of the 10 portfolios formed on momentum but with two different weighted methods, we have found that the mean-variance optimal set shifts to the left when we change from the average equal weighted method to the average value weighted method, but with the intersection near the expected monthly return of 0.018 and its standard deviation of 0.060. To increase the expected monthly return and therefore the associated risk of an efficient portfolio, the optimal portfolio weights for portfolios 2, 4, and 10 have to be increased for both weighted methods, in addtion to increasing the optimal portfolio weights for portfolios 6 and 9 for average equal weighted method.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\appendix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Residual Plots}

\newpage

\section{(AEW) Portfolio 1}

<<echo=FALSE>>=
portfolio <- 1
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 1}

<<echo=FALSE>>=
portfolio <- 1
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 2}

<<echo=FALSE>>=
portfolio <- 2
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 2}

<<echo=FALSE>>=
portfolio <- 2
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 3}

<<echo=FALSE>>=
portfolio <- 3
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 3}

<<echo=FALSE>>=
portfolio <- 3
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 4}

<<echo=FALSE>>=
portfolio <- 4
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 4}

<<echo=FALSE>>=
portfolio <- 4
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 5}

<<echo=FALSE>>=
portfolio <- 5
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 5}

<<echo=FALSE>>=
portfolio <- 5
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 6}

<<echo=FALSE>>=
portfolio <- 6
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 6}

<<echo=FALSE>>=
portfolio <- 6
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 7}

<<echo=FALSE>>=
portfolio <- 7
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 7}

<<echo=FALSE>>=
portfolio <- 7
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 8}

<<echo=FALSE>>=
portfolio <- 8
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 8}

<<echo=FALSE>>=
portfolio <- 8
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 9}

<<echo=FALSE>>=
portfolio <- 9
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 9}

<<echo=FALSE>>=
portfolio <- 9
residplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 10}

<<echo=FALSE>>=
portfolio <- 10
residplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 10}

<<echo=FALSE>>=
portfolio <- 10
residplotFun(data.l, 2, portfolio)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Plot of Cook's Distance and Covariance ratio}

\newpage

\section{Portfolio 1}

<<echo=FALSE>>=
toplotFun(data.l, 1, TRUE)
@

\newpage

\section{Portfolio 2}

<<echo=FALSE>>=
toplotFun(data.l, 2, TRUE)
@

\newpage

\section{Portfolio 3}

<<echo=FALSE>>=
toplotFun(data.l, 3, TRUE)
@

\newpage

\section{Portfolio 4}

<<echo=FALSE>>=
toplotFun(data.l, 4, TRUE)
@

\newpage

\section{Portfolio 5}

<<echo=FALSE>>=
toplotFun(data.l, 5, TRUE)
@

\newpage

\section{Portfolio 6}

<<echo=FALSE>>=
toplotFun(data.l, 6, TRUE)
@

\newpage

\section{Portfolio 7}

<<echo=FALSE>>=
toplotFun(data.l, 7, TRUE)
@

\newpage

\section{Portfolio 8}

<<echo=FALSE>>=
toplotFun(data.l, 8, TRUE)
@

\newpage

\section{Portfolio 9}

<<echo=FALSE>>=
toplotFun(data.l, 9, TRUE)
@

\newpage

\section{Portfolio 10}

<<echo=FALSE>>=
toplotFun(data.l, 10, TRUE)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Histogram of Cook's Distance and Covariance ratio}

\newpage

\section{Portfolio 1}

<<echo=FALSE>>=
toplotFun(data.l, 1, FALSE)
@

\newpage

\section{Portfolio 2}

<<echo=FALSE>>=
toplotFun(data.l, 2, FALSE)
@

\newpage

\section{Portfolio 3}

<<echo=FALSE>>=
toplotFun(data.l, 3, FALSE)
@

\newpage

\section{Portfolio 4}

<<echo=FALSE>>=
toplotFun(data.l, 4, FALSE)
@

\newpage

\section{Portfolio 5}

<<echo=FALSE>>=
toplotFun(data.l, 5, FALSE)
@

\newpage

\section{Portfolio 6}

<<echo=FALSE>>=
toplotFun(data.l, 6, FALSE)
@

\newpage

\section{Portfolio 7}

<<echo=FALSE>>=
toplotFun(data.l, 7, FALSE)
@

\newpage

\section{Portfolio 8}

<<echo=FALSE>>=
toplotFun(data.l, 8, FALSE)
@

\newpage

\section{Portfolio 9}

<<echo=FALSE>>=
toplotFun(data.l, 9, FALSE)
@

\newpage

\section{Portfolio 10}

<<echo=FALSE>>=
toplotFun(data.l, 10, FALSE)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Influence Measures of Influential Observations}

\newpage

\section{(AEW) Portfolio 1}

<<echo=FALSE>>=
portfolio <- 1
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 1}

<<echo=FALSE>>=
portfolio <- 1
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 2}

<<echo=FALSE>>=
portfolio <- 2
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 2}

<<echo=FALSE>>=
portfolio <- 2
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 3}

<<echo=FALSE>>=
portfolio <- 3
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 3}

<<echo=FALSE>>=
portfolio <- 3
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 4}

<<echo=FALSE>>=
portfolio <- 4
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 4}

<<echo=FALSE>>=
portfolio <- 4
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 5}

<<echo=FALSE>>=
portfolio <- 5
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 5}

<<echo=FALSE>>=
portfolio <- 5
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 6}

<<echo=FALSE>>=
portfolio <- 6
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 6}

<<echo=FALSE>>=
portfolio <- 6
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 7}

<<echo=FALSE>>=
portfolio <- 7
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 7}

<<echo=FALSE>>=
portfolio <- 7
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 8}

<<echo=FALSE>>=
portfolio <- 8
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 8}

<<echo=FALSE>>=
portfolio <- 8
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 9}

<<echo=FALSE>>=
portfolio <- 9
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 9}

<<echo=FALSE>>=
portfolio <- 9
exclusioninfoFun(data.v, portfolio) 
@

\newpage

\section{(AEW) Portfolio 10}

<<echo=FALSE>>=
portfolio <- 10
exclusioninfoFun(data.e, portfolio) 
@

\newpage

\section{(AVW) Portfolio 10}

<<echo=FALSE>>=
portfolio <- 10
exclusioninfoFun(data.v, portfolio) 
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Excluded Months By Panel and Weighting Method}

<<echo=FALSE>>=
list(
"(AEW) Panel B's Excluded Months" = exclusion.Be, 
"(AEW) Panel B's Number of Excluded Months" = length(exclusion.Be), 
"(AVW) Panel B's Excluded Months" = exclusion.Bv, 
"(AVW) Panel B's Number of Excluded Months" = length(exclusion.Bv), 
"(AEW) Panel C's Excluded Months" = 
 matrix(c(as.character(exclusion.Ce), rep("",3)), ncol = 5, byrow = TRUE), 
"(AEW) Panel C's Number of Excluded Months" = length(exclusion.Ce), 
"(AVW) Panel C's Excluded Months" = 
 matrix(c(as.character(exclusion.Cv), rep("",3)), ncol = 5, byrow = TRUE),
"(AVW) Panel C's Number of Excluded Months" = length(exclusion.Cv)
)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Plot of MOM vs Monthly Portfolio Return}

\newpage

\section{(AEW) Portfolio 1}

<<echo=FALSE>>=
portfolio <- 1
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 1}

<<echo=FALSE>>=
portfolio <- 1
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 2}

<<echo=FALSE>>=
portfolio <- 2
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 2}

<<echo=FALSE>>=
portfolio <- 2
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 3}

<<echo=FALSE>>=
portfolio <- 3
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 3}

<<echo=FALSE>>=
portfolio <- 3
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 4}

<<echo=FALSE>>=
portfolio <- 4
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 4}

<<echo=FALSE>>=
portfolio <- 4
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 5}

<<echo=FALSE>>=
portfolio <- 5
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 5}

<<echo=FALSE>>=
portfolio <- 5
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 6}

<<echo=FALSE>>=
portfolio <- 6
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 6}

<<echo=FALSE>>=
portfolio <- 6
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 7}

<<echo=FALSE>>=
portfolio <- 7
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 7}

<<echo=FALSE>>=
portfolio <- 7
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 8}

<<echo=FALSE>>=
portfolio <- 8
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 8}

<<echo=FALSE>>=
portfolio <- 8
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 9}

<<echo=FALSE>>=
portfolio <- 9
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 9}

<<echo=FALSE>>=
portfolio <- 9
momreturnplotFun(data.l, 2, portfolio)
@

\newpage

\section{(AEW) Portfolio 10}

<<echo=FALSE>>=
portfolio <- 10
momreturnplotFun(data.l, 1, portfolio)
@

\newpage

\section{(AVW) Portfolio 10}

<<echo=FALSE>>=
portfolio <- 10
momreturnplotFun(data.l, 2, portfolio)
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{R Codes}

<<howto, tidy=FALSE, eval=FALSE>>=
#############################################################################
### The Beginning of R codes for functions and variables
#############################################################################
### Import the data ###
#############################################################################

library(sandwich)
library(lmtest)
library(tseries)
library(xtable)
library(nortest)
library(car)

### Import the data
portfolios.e <- read.table("10_Portfolios_AveEquWeight.txt", header=T,
quote="\"")
portfolios.v <- read.table("10_Portfolios_AveValWeight.txt", header=T,
quote="\"")
researchfactor <- read.table("F-F_Research_Data_Factors.txt", header=T,
quote="\"")
momfactor <- read.table("F-F_Momentum_Factor.TXT", header=T,
quote="\"")

### Create function to combine the data
dataFun <- function(portfolios, researchfactor, momfactor){
    
    ### Combine the data
    data <- cbind(portfolios, researchfactor[,2:4], momfactor[,2],
    researchfactor[,5])
    
    ### Rename the variables
    names(data)[2:11] <- paste("Portfolio", 1:10, sep="")
    names(data)[12] <- "RMRF"
    names(data)[15:16] <- c("MOM", "RF")
    
    ### Convert the data in percentage
    data[,2:16] <- data[,2:16] / 100
    
    ### Compute the excess portfolio returns
    data[,2:11] <- data[,2:11] - data[,16]
    
    ### Convert the month from numeric to date format
    month.temp <- as.character(data[,1])
    data[,1] <- paste(substring(month.temp, 1, 4),
    substring(month.temp, 5, 6), "01", sep="-")
    data[,1] <- as.Date(data[,1])
    
    return (data)
}

### Get the average equal weighted (AEW) data with excess return
data.e <- dataFun(portfolios.e, researchfactor, momfactor)

### Get the average value weighted (AVW) data with excess return
data.v <- dataFun(portfolios.v, researchfactor, momfactor)

### Data with real return
pdata.e <- data.e
pdata.v <- data.v
pdata.e[,2:11] <- pdata.e[,2:11] + pdata.e[,16]
pdata.v[,2:11] <- pdata.v[,2:11] + pdata.v[,16]

### List of data with excess return
data.l <- list(data.e, data.v)

### List of data with real return
pdata.l <- list(pdata.e, pdata.v)

#############################################################################
### Create functions and variables for Chapter 1
#############################################################################

### Create function to generate residual plots
residplotFun <- function(datalist, datalistnum, portfolio){
  data <- datalist[[datalistnum]]
  i <- portfolio + 1
  fit <- lm(data[,i] ~ RMRF + SMB + HML + MOM, data = data)
  
  titlename <- paste(c("(AEW)", "(AVW)")[datalistnum],
                     "Residual Plots: Portfolio", portfolio, sep = " ")
  
  par(mfrow=c(3,2), oma=c(0,0,4,0), mar=c(5.1,4.1,1.1,2.1), las=1, cex=0.7)
  
  plot(data[,"Month"], residuals(fit), 
       main = "Time Series of Residuals",
       xlab = "Month", ylab = "Residuals",
       type = "l", col="blue")
  
  plot(fitted(fit), residuals(fit),
       main = "Fitted vs Residuals",
       xlab = "Fitted", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"RMRF"], residuals(fit),
       main = "RMRF vs Residuals",
       xlab = "RMRF", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"SMB"], residuals(fit),
       main = "SMB vs Residuals",
       xlab = "SMB", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"HML"], residuals(fit),
       main = "HML vs Residuals",
       xlab = "HML", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  plot(data[,"MOM"], residuals(fit),
       main = "MOM vs Residuals",
       xlab = "MOM", ylab = "Residuals", col = "blue")
  abline (h=0, col = "red")
  
  mtext(titlename, 3, 1, outer=TRUE, cex=1)
  
  # Restore default clipping rect
  par(mar = c(5, 4, 4, 2) + 0.1)
  par(mfrow=c(1,1))
}

#############################################################################

### Create function to compute p-value of Lilliefors test for normality
normtestpvalFun <- function(datalist){
  portfolio <- 1:10
  pvalAEW <- c(NA, 10)
  pvalAVW <- c(NA, 10)
  nullAEW <- c(NA, 10)
  nullAVW <- c(NA, 10)
  
  for(i in 1:10){
    fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[1]])
    resid1 <- residuals(fit1)
    pvalAEW[i] <- lillie.test(resid1)$p.value
    if(pvalAEW[i] < 0.05){
      nullAEW[i] <- "Reject"
    } else {nullAEW[i] <- "Accept"}
    
    fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[2]])
    resid2 <- residuals(fit2)
    pvalAVW[i] <- lillie.test(resid2)$p.value
    if(pvalAVW[i] < 0.05){
      nullAVW[i] <- "Reject"
    } else {nullAVW[i] <- "Accept"}
  }
  
  return(data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW))
}

### Compute the p-value
normtestpval.A <- normtestpvalFun(data.l)

#############################################################################

### Create function to compute mean of residuals
meanresidFun <- function(datalist){
  portfolio <- 1:10
  meanresidAEW <- c(NA, 10)
  meanresidAVW <- c(NA, 10)
  
  for(i in 1:10){
    fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[1]])
    meanresidAEW[i] <- mean(residuals(fit1))
    
    fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, data=datalist[[2]])
    meanresidAVW[i] <- mean(residuals(fit2))
  }
  
  return(data.frame(portfolio, meanresidAEW, meanresidAVW))
}

### Compute the mean of residuals
meanresid.A <- meanresidFun(data.l)

#############################################################################

### Compute the p-values of non-constant variance score test
datalist <- data.l
portfolio <- 1:10
pvalAEW <- c(NA, 10)
pvalAVW <- c(NA, 10)
nullAEW <- c(NA, 10)
nullAVW <- c(NA, 10)

for(i in 1:10){
  fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[1]])
  pvalAEW[i] <- ncvTest(fit1)$p
  if(pvalAEW[i] < 0.05){
    nullAEW[i] <- "Reject"
  } else {nullAEW[i] <- "Accept"}
  
  fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[2]])
  pvalAVW[i] <- ncvTest(fit2)$p
  if(pvalAVW[i] < 0.05){
    nullAVW[i] <- "Reject"
  } else {nullAVW[i] <- "Accept"}
}

ncvtestpval.A <- data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW)
rm(datalist, portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW, i, fit1, fit2)

#############################################################################
### Create functions and variables for Chapter 2
#############################################################################

### Create function to plot Cook's distance
CiPlotFun <- function(datalist, datalistnum, portfolio, toplot = TRUE){
    data <- datalist[[datalistnum]]
    colnum <- portfolio + 1
    
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    
    ti <- as.matrix(rep(1,N))
    Ci <- as.matrix(rep(1,N))
    for(j in 1:N){
        ti[j] <- resid[j] / (sqrt(sigma2hat) * sqrt(1 - P[j,j]))
        Ci[j] <- ( ((ti[j])^2) / p ) * ( P[j,j] / (1 - P[j,j]) )
    }
    
    calpt.find <- Ci > 1
    Month <- as.matrix(data[,1])
    Cidata <- cbind(Month, Ci)
    colnames(Cidata) <- c("Month", "Ci")
    Cidata <- as.data.frame(Cidata)
    Cidata[,1] <- as.Date(Cidata[,1])
    calpt.count <- sum(calpt.find)
    Cidata.sort <- Cidata[order(-Ci),]
    
    dataname <- c("(AEW)", "(AVW)")[datalistnum]
    if(toplot == TRUE){
        plot(Cidata[,"Month"], Cidata[,"Ci"],
        main = paste(dataname, "\nCook's Distance: Portfolio", portfolio,
        sep = " "),
        xlab = "Month", ylab = "Cook's Distance", col="blue")
        points(Cidata[,"Month"][calpt.find], Cidata[,"Ci"][calpt.find],
        col="red", pch=19)
    } else {
        hist(Cidata[,"Ci"], breaks = 100,
        col = "skyblue", xlab = "Cook's Distance",
        main = paste(dataname, "Histogram\nCook's Distance: Portfolio",
        portfolio, sep = " "))
    }
}

#############################################################################

### Create function to plot covariance ratio
CVRiPlotFun <- function(datalist, datalistnum, portfolio, toplot = TRUE){
    data <- datalist[[datalistnum]]
    colnum <- portfolio + 1
    
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    sigma2hat.i <- as.matrix(rep(1,N))
    for(j in 1:N){
        sigma2hat.i[j] <- (N - p)*sigma2hat/(N - p -1) -
        resid[j]^2/((N - p - 1)*(1 - P[j,j]))
    }
    
    CVRi <- as.matrix(rep(1,N))
    for(j in 1:N){
        CVRi[j] <- ((sigma2hat.i[j]/sigma2hat)^p) / (1 - P[j,j])
    }
    
    calpt.find <- abs(CVRi - 1 ) > 3*p/N
    abs.CVRi_1 <- abs(CVRi - 1 )
    calpt.crit <- 3*p/N
    Month <- as.matrix(data[,1])
    CVRidata <- data.frame(Month, CVRi, abs.CVRi_1, calpt.crit)
    CVRidata[,1] <- as.Date(CVRidata[,1])
    
    calpt.count <- sum(calpt.find)
    CVRidata.sort <- CVRidata[order(-abs.CVRi_1),]
    
    dataname <- c("(AEW)", "(AVW)")[datalistnum]
    if(toplot == TRUE){
        plot(CVRidata$Month, CVRidata$CVRi,
        main = paste(dataname, "\nCovariance Ratio: Portfolio", portfolio,
        sep = " "),
        xlab = "Month", ylab = "Covariance Ratio", col="blue")
        points(CVRidata$Month[calpt.find], CVRidata$CVRi[calpt.find], col="red",
        pch=19)
    } else {
        hist(CVRidata$CVRi, breaks = 100,
        col = "skyblue", xlab = "Covariance Ratio",
        main = paste(dataname, "Histogram\nCovariance Ratio: Portfolio",
        portfolio, sep = " "))
    }
}

#############################################################################

### Create function to combine plots of Cook's distance and covariance ratio
toplotFun <- function(datalist, portfolio, toplot = TRUE){
    if(toplot == TRUE){
        par(mfcol=c(2,2))
        for(i in 1:2){
            CiPlotFun(datalist, i, portfolio)
            CVRiPlotFun(datalist, i, portfolio)
        }
        par(mfrow=c(1,1))
    } else {
        par(mfcol=c(4,1))
        for(i in 1:2){
            CiPlotFun(datalist, i, portfolio, FALSE)
            
        }
        for(i in 1:2){
            CVRiPlotFun(datalist, i, portfolio, FALSE)
        }
        par(mfrow=c(1,1))
    }
}

#############################################################################

### Create function to plot MOM vs monthly portfolio returns
momreturnplotFun <- function(datalist, datalistnum, portfolio){
    data <- datalist[[datalistnum]]
    i <- c(1:7)
    i <- 130 * i
    start <- c(1, i + 1)
    end <- c(i, nrow(data))
    
    j <- portfolio + 1
    temp <- data[,j] + data[,16]
    ymin <- min(min(temp), min(data[,15]))
    ymax <- max(max(temp), max(data[,15]))
    
    titlename <- paste(c("(AEW)", "(AVW)")[datalistnum],
    "MOM & Monthly Portfolio Return: Portfolio", portfolio,
    sep = " ")
    
    par(mfrow=c(4,2), oma=c(3,3,4,0), mar=c(4,2,1,1), las=1, cex=0.7)
    for(k in 1:8){
        plot(data[start[k]:end[k],1], temp[start[k]:end[k]],
        xlab = "", ylab = "", ylim=c(ymin,ymax), type="l", col = "blue")
        lines(data[start[k]:end[k],1], data[start[k]:end[k],15], type="l",
        col = "red")
        legend("topleft", c("Return", "MOM"), cex=0.5, col=c("blue", "red"),
        lty=1, bty="n")
    }
    mtext("Month", 1, 1, outer=TRUE, cex = 0.75)
    mtext("MOM & Monthly Portfolio Return", 2, 1, outer=TRUE, las=0,
    cex = 0.75)
    mtext(titlename, 3, 1, outer=TRUE, cex=1)
    
    # Restore default clipping rect
    par(mar = c(5, 4, 4, 2) + 0.1)
    par(mfrow=c(1,1))
}

#############################################################################

### Create function to compute the Cook's distance for Influential Points
CiCalptFun <- function(data, colnum){
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    
    ti <- as.matrix(rep(1,N))
    Ci <- as.matrix(rep(1,N))
    for(j in 1:N){
        ti[j] <- resid[j] / (sqrt(sigma2hat) * sqrt(1 - P[j,j]))
        Ci[j] <- ( ((ti[j])^2) / p ) * ( P[j,j] / (1 - P[j,j]) )
    }
    
    calpt.find <- Ci > 1
    Month <- as.matrix(data[,1])
    Cidata <- cbind(Month, Ci)
    colnames(Cidata) <- c("Month", "Ci")
    Cidata <- as.data.frame(Cidata)
    Cidata[,1] <- as.Date(Cidata[,1])
    calpt.count <- sum(calpt.find)
    Cidata.sort <- Cidata[order(-Ci),]
    
    if(calpt.count == 0){
        return(NULL)
    } else {
        return(Cidata.sort[1:calpt.count,])
    }
    
}

#############################################################################

### Create function to compute the covariance ratios for Influential Points
CVRiCalptFun <- function(data, colnum){
    N <- nrow(data)
    X <- as.matrix(data[,12:15])
    X <- cbind(as.matrix(rep(1,N)), X)
    p <- ncol(X)
    Y <- as.matrix(data[,colnum, drop=FALSE])
    betahat <- solve(t(X) %*% X) %*% t(X) %*% Y
    Yhat <- X %*% betahat
    P <- X %*% solve(t(X) %*% X) %*% t(X)
    resid <- Y - Yhat
    sigma2hat <- (t(resid) %*% resid) / (N - p)
    sigma2hat.i <- as.matrix(rep(1,N))
    
    for(j in 1:N){
        sigma2hat.i[j] <- (N - p)*sigma2hat/(N - p -1) -
        resid[j]^2/((N - p - 1)*(1 - P[j,j]))
    }
    
    CVRi <- as.matrix(rep(1,N))
    for(j in 1:N){
        CVRi[j] <- ((sigma2hat.i[j]/sigma2hat)^p) / (1 - P[j,j])
    }
    
    calpt.find <- abs(CVRi - 1 ) > 3*p/N
    abs.CVRi_1 <- abs(CVRi - 1 )
    calpt.crit <- 3*p/N
    Month <- as.matrix(data[,1])
    CVRidata <- data.frame(Month, CVRi, abs.CVRi_1, calpt.crit)
    CVRidata[,1] <- as.Date(CVRidata[,1])
    
    calpt.count <- sum(calpt.find)
    CVRidata.sort <- CVRidata[order(-abs.CVRi_1),]
    
    if(calpt.count == 0){
        return(NULL)
    } else {
        return(CVRidata.sort[1:calpt.count,])
    }
}

#############################################################################

### Create function to compute influence measures of influential observations
exclusioninfoFun <- function(data, portfolio){
    Ci <- as.list(c(1:10))
    CVARi <- as.list(c(1:10))
    
    for(i in 1:10){
        
        if(is.null(CiCalptFun(data, i+1)) == TRUE){
            Ci[[i]] <- "None"
        } else {
            Ci[[i]] <- CiCalptFun(data, i+1)
        }
        
        if(is.null(CVRiCalptFun(data, i+1)) == TRUE){
            CVARi[[i]] <- "None"
        } else {
            CVARi[[i]] <- CVRiCalptFun(data, i+1)
        }
        
    }
    
    exclusion.info <- list("Cook's Distance" = Ci[[portfolio]],
    "Covariance Ratio" = CVARi[[portfolio]])
    
    return(exclusion.info)
}

#############################################################################
### Create functions and variables for Chapter 3
#############################################################################

### Get the excluded months
exclusion.Be <- c()
exclusion.Bv <- c()
exclusion.Ce <- c()
exclusion.Cv <- c()

for(i in 1:10){
    exclusion.Be <- c(exclusion.Be,
    CiCalptFun(data.e, i+1)[,1])
    exclusion.Bv <- c(exclusion.Bv,
    CiCalptFun(data.v, i+1)[,1])
    exclusion.Ce <- c(exclusion.Ce,
    CVRiCalptFun(data.e, i+1)[,1])
    exclusion.Cv <- c(exclusion.Cv,
    CVRiCalptFun(data.v, i+1)[,1])
}

exclusion.Be <- as.Date(sort(unique(exclusion.Be)))
exclusion.Bv <- as.Date(sort(unique(exclusion.Bv)))
exclusion.Ce <- as.Date(sort(unique(exclusion.Ce)))
exclusion.Cv <- as.Date(sort(unique(exclusion.Cv)))

### Get the subset lists
subset.Ae <- as.list(c(1:10))
subset.Av <- as.list(c(1:10))
subset.Be <- as.list(c(1:10))
subset.Bv <- as.list(c(1:10))
subset.Ce <- as.list(c(1:10))
subset.Cv <- as.list(c(1:10))

for(i in 1:10){
    subset.Ae[[i]] <- rep(TRUE, nrow(data.e))
    subset.Av[[i]] <- rep(TRUE, nrow(data.v))
    subset.Be[[i]] <- !(data.e[,1] %in% exclusion.Be)
    subset.Bv[[i]] <- !(data.v[,1] %in% exclusion.Bv)
    subset.Ce[[i]] <- !(data.e[,1] %in% exclusion.Ce)
    subset.Cv[[i]] <- !(data.v[,1] %in% exclusion.Cv)
}

### Panel B data
  data.Be <- data.e[subset.Be[[1]],] ### AEW excess return
  data.Bv <- data.v[subset.Bv[[1]],] ### AVW excess return
  pdata.Be <- pdata.e[subset.Be[[1]],] ### AEW real return
  pdata.Bv <- pdata.v[subset.Bv[[1]],] ### AVW real return
  
  ### List of data
  data.Bl <- list(data.Be, data.Bv)
  pdata.Bl <- list(pdata.Be, pdata.Bv)

### Panel C data 
  data.Ce <- data.e[subset.Ce[[1]],] ### AEW excess return
  data.Cv <- data.v[subset.Cv[[1]],] ### AVW excess return
  pdata.Ce <- pdata.e[subset.Ce[[1]],] ### AEW real return
  pdata.Cv <- pdata.v[subset.Cv[[1]],] ### AVW real return
  
  ### List of data
  data.Cl <- list(data.Ce, data.Cv)
  pdata.Cl <- list(pdata.Ce, pdata.Cv)

#############################################################################

### Create function to compute the statistics
olsFun <- function(data, portfolio){
  varindex <- portfolio + 1
  fit <- lm(data[,varindex] ~ RMRF + SMB + HML + MOM, data=data)
  
  varnames <- c("Intercept", "RMRF", "SMB", "HML", "MOM")
  beta <- matrix(NA, 1, 5)
  beta[1,] <- fit$coefficients
  colnames(beta) <- varnames
  
  beta.t <- matrix(NA, 1, 5)
  beta.t[1,] <- coeftest(fit)[,"t value"]
  colnames(beta.t) <- varnames
  
  beta.NWt <- matrix(NA, 1, 5)
  beta.NWt[1,] <- coeftest(fit, NeweyWest(fit))[,"t value"]
  colnames(beta.NWt) <- varnames
  
  return(list("Beta" = beta, "Beta t-statistic" = beta.t,
              "Beta Newey-West t-statistic" = beta.NWt))
}

#############################################################################

### Create function to generate the panel summary
panelsummaryFun <- function(data, subset){
    panel.sum <- olsFun(data[subset[[1]], ], 1)
    for(i in 2:10){
        stats.temp <- olsFun(data[subset[[i]], ], i)
        for(j in 1:3){
            panel.sum[[j]] <- rbind(panel.sum[[j]], stats.temp[[j]])
        }
    }
    
    Portfolio <- 1:10
    for(j in 1:3){
        panel.sum[[j]] <- cbind(Portfolio, panel.sum[[j]])
    }
    
    return(panel.sum)
}

#############################################################################

### Generate summary for each panel
### AEW
panelAe.summary <- panelsummaryFun(data.e, subset.Ae)
panelBe.summary <- panelsummaryFun(data.e, subset.Be)
panelCe.summary <- panelsummaryFun(data.e, subset.Ce)

### AVW
panelAv.summary <- panelsummaryFun(data.v, subset.Av)
panelBv.summary <- panelsummaryFun(data.v, subset.Bv)
panelCv.summary <- panelsummaryFun(data.v, subset.Cv)

#############################################################################

### (AEW) Beta
### Panel A
beta.Ae <- panelAe.summary[[1]]

### Panel B
beta.Be <- panelBe.summary[[1]]

### Percentage change from panel A to panel B
beta.ABe <- ((panelBe.summary[[1]] - panelAe.summary[[1]]) /
panelAe.summary[[1]]) * 100
beta.ABe[,1] <- 1:10

#############################################################################

### (AVW) Beta
### Panel A
beta.Av <- panelAv.summary[[1]]

### Panel B
beta.Bv <- panelBv.summary[[1]]

### Percentage change from panel A to panel B
beta.ABv <- ((panelBv.summary[[1]] - panelAv.summary[[1]]) /
panelAv.summary[[1]]) * 100
beta.ABv[,1] <- 1:10

#############################################################################

### (AEW) Beta t-statistic
### Panel A
beta_t.Ae <- panelAe.summary[[2]]

### Panel C
beta_t.Ce <- panelCe.summary[[2]]

### Percentage change from panel A to panel C
beta_t.ACe <- ((panelCe.summary[[2]] - panelAe.summary[[2]]) /
panelAe.summary[[2]]) * 100
beta_t.ACe[,1] <- 1:10

#############################################################################

### (AEW) Beta Newey-West t-statistic
### Panel A
beta_nwt.Ae <- panelAe.summary[[3]]

### Panel C
beta_nwt.Ce <- panelCe.summary[[3]]

### Percentage change from panel A to panel C
beta_nwt.ACe <- ((panelCe.summary[[3]] - panelAe.summary[[3]]) /
panelAe.summary[[3]]) * 100
beta_nwt.ACe[,1] <- 1:10

#############################################################################

### (AVW) Beta t-statistic
### Panel A
beta_t.Av <- panelAv.summary[[2]]

### Panel C
beta_t.Cv <- panelCv.summary[[2]]

### Percentage change from panel A to panel C
beta_t.ACv <- ((panelCv.summary[[2]] - panelAv.summary[[2]]) /
panelAv.summary[[2]]) * 100
beta_t.ACv[,1] <- 1:10

#############################################################################

### (AVW) Beta Newey-West t-statistic
### Panel A
beta_nwt.Av <- panelAv.summary[[3]]

### Panel C
beta_nwt.Cv <- panelCv.summary[[3]]

### Percentage change from panel A to panel C
beta_nwt.ACv <- ((panelCv.summary[[3]] - panelAv.summary[[3]]) /
panelAv.summary[[3]]) * 100
beta_nwt.ACv[,1] <- 1:10

#############################################################################

### Examine the OLS assumptions for panels B and C

### Panel B

### Compute the p-Values of Lilliefors test for normality
normtestpval.B <- normtestpvalFun(data.Bl)

### Compute the p-values of non-constant variance score test
datalist <- data.Bl
portfolio <- 1:10
pvalAEW <- c(NA, 10)
pvalAVW <- c(NA, 10)
nullAEW <- c(NA, 10)
nullAVW <- c(NA, 10)

for(i in 1:10){
  fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[1]])
  pvalAEW[i] <- ncvTest(fit1)$p
  if(pvalAEW[i] < 0.05){
    nullAEW[i] <- "Reject"
  } else {nullAEW[i] <- "Accept"}
  
  fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[2]])
  pvalAVW[i] <- ncvTest(fit2)$p
  if(pvalAVW[i] < 0.05){
    nullAVW[i] <- "Reject"
  } else {nullAVW[i] <- "Accept"}
}

ncvtestpval.B <- data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW)
rm(datalist, portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW, i, fit1, fit2)

### Panel C

### Compute the p-Values of Lilliefors test for normality
normtestpval.C <-normtestpvalFun(data.Cl)

### Compute the p-values of non-constant variance score test
datalist <- data.Cl
portfolio <- 1:10
pvalAEW <- c(NA, 10)
pvalAVW <- c(NA, 10)
nullAEW <- c(NA, 10)
nullAVW <- c(NA, 10)

for(i in 1:10){
  fit1 <- lm(datalist[[1]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[1]])
  pvalAEW[i] <- ncvTest(fit1)$p
  if(pvalAEW[i] < 0.05){
    nullAEW[i] <- "Reject"
  } else {nullAEW[i] <- "Accept"}
  
  fit2 <- lm(datalist[[2]][,i+1] ~ RMRF + SMB + HML + MOM, 
             data=datalist[[2]])
  pvalAVW[i] <- ncvTest(fit2)$p
  if(pvalAVW[i] < 0.05){
    nullAVW[i] <- "Reject"
  } else {nullAVW[i] <- "Accept"}
}

ncvtestpval.C <- data.frame(portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW)
rm(datalist, portfolio, pvalAEW, pvalAVW, nullAEW, nullAVW, i, fit1, fit2)


#############################################################################
### Create functions and variables for Chapter 4
#############################################################################

### Create function to compute mean and sd of monthly portfolio return

statFun <- function(pdata, subset){
    statsummary <- matrix(NA, 10, 3)
    rownames(statsummary) <- rep("", 10)
    colnames(statsummary) <- c("Portfolio", "Mean", "SD")
    statsummary[,1] <- 1:10
    for(i in 1:10){
        statsummary[i,2] <- mean(pdata[subset[[i]], i+1])
        statsummary[i,3] <- sd(pdata[subset[[i]], i+1])
    }
    return(statsummary)
}

#############################################################################

### (AEW) Mean and sd of return
### Panel A
meansd.Ae <- statFun(pdata.e, subset.Ae)

### Panel B
meansd.Be <- statFun(pdata.e, subset.Be)

### Panel C
meansd.Ce <- statFun(pdata.e, subset.Ce)

### Percentage change from panel A to panel B
meansd.ABe <- ((meansd.Be - meansd.Ae) / meansd.Ae) * 100
meansd.ABe[,1] <- 1:10

### Percentage change from panel A to panel C
meansd.ACe <- ((meansd.Ce - meansd.Ae) / meansd.Ae) * 100
meansd.ACe[,1] <- 1:10

#############################################################################

### (AVW) Mean and sd of return
### Panel A
meansd.Av <- statFun(pdata.v, subset.Av)

### Panel B
meansd.Bv <- statFun(pdata.v, subset.Bv)

### Panel C
meansd.Cv <- statFun(pdata.v, subset.Cv)

### Percentage change from panel A to panel B
meansd.ABv <- ((meansd.Bv - meansd.Av) / meansd.Av) * 100
meansd.ABv[,1] <- 1:10

### Percentage change from panel A to panel C
meansd.ACv <- ((meansd.Cv - meansd.Av) / meansd.Av) * 100
meansd.ACv[,1] <- 1:10

#############################################################################

### Create function to compute numbers for mean-variance optimal portfolio

effpfFun <- function(pdata, mu.min, mu.max, nobs, toshorts = TRUE){
    
    pfoptFun <- function(mu, pdata, toshorts){
        pfmu <- matrix(colMeans(pdata[,2:11]),nrow=1)
        pfcov <- cov(pdata[,2:11])
        if(toshorts == TRUE){
            pfopt.sol <- portfolio.optim(x = pfmu, pm = mu,
            covmat = pfcov, shorts = toshorts)
        } else {
            pfopt.sol <- portfolio.optim(x = pfmu, pm = mu,
            covmat = pfcov, shorts = toshorts,
            reslow = rep(0,10), reshigh = rep(1,10))
        }
        
        wt <- pfopt.sol$pw
        SD <- sqrt((t(wt) %*% pfcov %*% wt)[1])
        return(c(mu, SD, wt))
    }
    
    mu.list <- seq(mu.min, mu.max, length.out = nobs)
    effpf <- t(sapply(mu.list, pfoptFun, pdata, toshorts))
    colnames(effpf) <- c("Mean", "SD", paste("W", 1:10, sep=""))
    rownames(effpf) <- rep("", nrow(effpf))
    
    return(effpf)
}

#############################################################################

### (AEW) Mean-Variance Optimal Portfolio

### Mean of each monthly portfolio return
pfmu.e <- matrix(colMeans(pdata.e[,2:11]),nrow=1)

### Max and min of the mean
pmin.e <- min(pfmu.e) + 1e-4
pmax.e <- max(pfmu.e) - 1e-4

### Mean-Variance Optimal Portfolio
effpfT.e <- effpfFun(pdata.e, 0, 0.5, 1000, TRUE)
effpfF.e <- effpfFun(pdata.e, pmin.e, pmax.e, 100, FALSE)

### Subset for nicer layout
min_x.e <- min(effpfT.e[,2], effpfF.e[,2])
max_x.e <- max(effpfF.e[,2])
effpfT.e <- effpfT.e[effpfT.e[,2] <= max_x.e,]

#############################################################################

### (AVW) Mean-Variance Optimal Portfolio

### Mean of each monthly portfolio return
pfmu.v <- matrix(colMeans(pdata.v[,2:11]),nrow=1)

### Max and min of the mean
pmin.v <- min(pfmu.v) + 1e-4
pmax.v <- max(pfmu.v) - 1e-4

### Mean-Variance Optimal Portfolio
effpfT.v <- effpfFun(pdata.v, 0, 0.5, 1000, TRUE)
effpfF.v <- effpfFun(pdata.v, pmin.v, pmax.v, 100, FALSE)

### Subset for nicer layout
min_x.v <- min(effpfT.v[,2], effpfF.v[,2])
max_x.v <- max(effpfF.v[,2])
effpfT.v <- effpfT.v[effpfT.v[,2] <= max_x.v,]

#############################################################################

### (Shorts Allowed) Mean-Variance Optimal Portfolio: AEW vs AVW
pmin <- 0
pmax <- 0.03

effpf.e <- effpfFun(pdata.e, pmin, pmax, 31, TRUE)
effpf.v <- effpfFun(pdata.v, pmin, pmax, 31, TRUE)

min_x.ev <- min(effpf.e[,2], effpf.v[,2])
max_x.ev <- max(effpf.e[,2], effpf.v[,2])

#############################################################################
### The End of R codes for functions and variables
#############################################################################

#############################################################################
### The Beginning of R codes for Chapter
#############################################################################
### Chapter 1: Introduction
#############################################################################

### OLS Assumptions

### Residual Plots
### (e.g. AEW Portfolio 2)
residplotFun(data.l, 1, 2)

### P-Values of Lilliefors Test for Normality
normtestpval.A

### Mean of Residuals
meanresid.A

### P-values of Non-constant Variance Score Test
ncvtestpval.A

#############################################################################
### Chapter 2: Identification of Influential Observation
#############################################################################

### Results from Cook's Distance and Covariance Ratio

### Plot of Cook's Distance and Covariance Ratio
### (e.g. Portfolio 2)
toplotFun(data.l, 2, TRUE)

### Histogram of Cook's Distance and Covariance Ratio
### (e.g. Portfolio 2)
toplotFun(data.l, 2, FALSE)

### Plot of MOM vs Monthly Portfolio Return
### (e.g. AEW Portfolio 2)
momreturnplotFun(data.l, 1, 2)


#############################################################################
### Chapter 3: Influence of Influential Observation
#############################################################################

### (AEW) Plots of Beta

par(mfrow=c(3,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta.Ae[,2], beta.Be[,2])
max_y <- max(beta.Ae[,2], beta.Be[,2])
plot(beta.Ae[,1], beta.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Intercept",
main = "(AEW) Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Ae[,3], beta.Be[,3])
max_y <- max(beta.Ae[,3], beta.Be[,3])
plot(beta.Ae[,1], beta.Ae[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(RMRF)",
main = "(AEW) Beta of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Ae[,4], beta.Be[,4])
max_y <- max(beta.Ae[,4], beta.Be[,4])
plot(beta.Ae[,1], beta.Ae[,4], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(SMB)",
main = "(AEW) Beta of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Ae[,5], beta.Be[,5])
max_y <- max(beta.Ae[,5], beta.Be[,5])
plot(beta.Ae[,1], beta.Ae[,5], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(HML)",
main = "(AEW) Beta of HML")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Ae[,6], beta.Be[,6])
max_y <- max(beta.Ae[,6], beta.Be[,6])
plot(beta.Ae[,1], beta.Ae[,6], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(MOM)",
main = "(AEW) Beta of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta.Be[,1], beta.Be[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.ABe[,2:6])
max_y <- max(beta.ABe[,2:6])
plot(beta.ABe[,1], beta.ABe[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AEW) % Change from Panel A to Panel B")
axis(1, at=1:10, lab=1:10)
lines(beta.ABe[,1], beta.ABe[,3], type="o", col="green")
lines(beta.ABe[,1], beta.ABe[,4], type="o", col="blue")
lines(beta.ABe[,1], beta.ABe[,5], type="o", col="purple")
lines(beta.ABe[,1], beta.ABe[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y,
c("Intercept", "Beta(RMRF)", "Beta(SMB)", "Beta(HML)", "Beta(MOM)"),
cex=0.6, col=c("red", "green", "blue", "purple", "brown"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### (AVW) Plots of Beta

par(mfrow=c(3,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta.Av[,2], beta.Bv[,2])
max_y <- max(beta.Av[,2], beta.Bv[,2])
plot(beta.Av[,1], beta.Av[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Intercept",
main = "(AVW) Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Av[,3], beta.Bv[,3])
max_y <- max(beta.Av[,3], beta.Bv[,3])
plot(beta.Av[,1], beta.Av[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(RMRF)",
main = "(AVW) Beta of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Av[,4], beta.Bv[,4])
max_y <- max(beta.Av[,4], beta.Bv[,4])
plot(beta.Av[,1], beta.Av[,4], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(SMB)",
main = "(AVW) Beta of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Av[,5], beta.Bv[,5])
max_y <- max(beta.Av[,5], beta.Bv[,5])
plot(beta.Av[,1], beta.Av[,5], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(HML)",
main = "(AVW) Beta of HML")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.Av[,6], beta.Bv[,6])
max_y <- max(beta.Av[,6], beta.Bv[,6])
plot(beta.Av[,1], beta.Av[,6], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Beta(MOM)",
main = "(AVW) Beta of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta.Bv[,1], beta.Bv[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel B"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta.ABv[,2:6])
max_y <- max(beta.ABv[,2:6])
plot(beta.ABv[,1], beta.ABv[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AVW) % Change from Panel A to Panel B")
axis(1, at=1:10, lab=1:10)
lines(beta.ABv[,1], beta.ABv[,3], type="o", col="green")
lines(beta.ABv[,1], beta.ABv[,4], type="o", col="blue")
lines(beta.ABv[,1], beta.ABv[,5], type="o", col="purple")
lines(beta.ABv[,1], beta.ABv[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y,
c("Intercept", "Beta(RMRF)", "Beta(SMB)", "Beta(HML)", "Beta(MOM)"),
cex=0.6, col=c("red", "green", "blue", "purple", "brown"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of Intercept t-statistic

par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,2], beta_t.Ce[,2])
max_y <- max(beta_t.Ae[,2], beta_t.Ce[,2])
plot(beta_t.Ae[,1], beta_t.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(Intercept)",
main = "(AEW) T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_t.Av[,2], beta_t.Cv[,2])
max_y <- max(beta_t.Av[,2], beta_t.Cv[,2])
plot(beta_t.Av[,1], beta_t.Av[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(Intercept)",
main = "(AVW) T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,2], beta_nwt.Ce[,2])
max_y <- max(beta_nwt.Ae[,2], beta_nwt.Ce[,2])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(Intercept)",
main = "(AEW) NW T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Av[,2], beta_nwt.Cv[,2])
max_y <- max(beta_nwt.Av[,2], beta_nwt.Cv[,2])
plot(beta_nwt.Av[,1], beta_nwt.Av[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(Intercept)",
main = "(AVW) NW T-Stats of Intercept")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,2], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of RMRF t-statistic

par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,3], beta_t.Ce[,3])
max_y <- max(beta_t.Ae[,3], beta_t.Ce[,3])
plot(beta_t.Ae[,1], beta_t.Ae[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(RMRF)",
main = "(AEW) T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_t.Av[,3], beta_t.Cv[,3])
max_y <- max(beta_t.Av[,3], beta_t.Cv[,3])
plot(beta_t.Av[,1], beta_t.Av[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(RMRF)",
main = "(AVW) T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,3], beta_nwt.Ce[,3])
max_y <- max(beta_nwt.Ae[,3], beta_nwt.Ce[,3])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(RMRF)",
main = "(AEW) NW T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Av[,3], beta_nwt.Cv[,3])
max_y <- max(beta_nwt.Av[,3], beta_nwt.Cv[,3])
plot(beta_nwt.Av[,1], beta_nwt.Av[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(RMRF)",
main = "(AVW) NW T-Stats of RMRF")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,3], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of SMB t-statistic

par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,4], beta_t.Ce[,4])
max_y <- max(beta_t.Ae[,4], beta_t.Ce[,4])
plot(beta_t.Ae[,1], beta_t.Ae[,4], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(SMB)",
main = "(AEW) T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_t.Av[,4], beta_t.Cv[,4])
max_y <- max(beta_t.Av[,4], beta_t.Cv[,4])
plot(beta_t.Av[,1], beta_t.Av[,4], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(SMB)",
main = "(AVW) T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,4], beta_nwt.Ce[,4])
max_y <- max(beta_nwt.Ae[,4], beta_nwt.Ce[,4])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,4], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(SMB)",
main = "(AEW) NW T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Av[,4], beta_nwt.Cv[,4])
max_y <- max(beta_nwt.Av[,4], beta_nwt.Cv[,4])
plot(beta_nwt.Av[,1], beta_nwt.Av[,4], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(SMB)",
main = "(AVW) NW T-Stats of SMB")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,4], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of HML t-statistic

par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,5], beta_t.Ce[,5])
max_y <- max(beta_t.Ae[,5], beta_t.Ce[,5])
plot(beta_t.Ae[,1], beta_t.Ae[,5], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(HML)",
main = "(AEW) T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_t.Av[,5], beta_t.Cv[,5])
max_y <- max(beta_t.Av[,5], beta_t.Cv[,5])
plot(beta_t.Av[,1], beta_t.Av[,5], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(HML)",
main = "(AVW) T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,5], beta_nwt.Ce[,5])
max_y <- max(beta_nwt.Ae[,5], beta_nwt.Ce[,5])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,5], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(HML)",
main = "(AEW) NW T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Av[,5], beta_nwt.Cv[,5])
max_y <- max(beta_nwt.Av[,5], beta_nwt.Cv[,5])
plot(beta_nwt.Av[,1], beta_nwt.Av[,5], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(HML)",
main = "(AVW) NW T-Stats of HML")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,5], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of MOM t-statistic

par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.Ae[,6], beta_t.Ce[,6])
max_y <- max(beta_t.Ae[,6], beta_t.Ce[,6])
plot(beta_t.Ae[,1], beta_t.Ae[,6], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(MOM)",
main = "(AEW) T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Ce[,1], beta_t.Ce[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_t.Av[,6], beta_t.Cv[,6])
max_y <- max(beta_t.Av[,6], beta_t.Cv[,6])
plot(beta_t.Av[,1], beta_t.Av[,6], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "T(MOM)",
main = "(AVW) T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_t.Cv[,1], beta_t.Cv[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Ae[,6], beta_nwt.Ce[,6])
max_y <- max(beta_nwt.Ae[,6], beta_nwt.Ce[,6])
plot(beta_nwt.Ae[,1], beta_nwt.Ae[,6], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(MOM)",
main = "(AEW) NW T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Ce[,1], beta_nwt.Ce[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

min_y <- min(beta_nwt.Av[,6], beta_nwt.Cv[,6])
max_y <- max(beta_nwt.Av[,6], beta_nwt.Cv[,6])
plot(beta_nwt.Av[,1], beta_nwt.Av[,6], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "NWT(MOM)",
main = "(AVW) NW T-Stats of MOM")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.Cv[,1], beta_nwt.Cv[,6], type="o", col="red")
legend(10.5, max_y, c("Panel A", "Panel C"), cex=0.6, col=c("blue", "red"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of % change in t-statistic

par(mfrow=c(2,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,5))

min_y <- min(beta_t.ACe[,2:6])
max_y <- max(beta_t.ACe[,2:6])
plot(beta_t.ACe[,1], beta_t.ACe[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AEW) % Change in T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_t.ACe[,1], beta_t.ACe[,3], type="o", col="green")
lines(beta_t.ACe[,1], beta_t.ACe[,4], type="o", col="blue")
lines(beta_t.ACe[,1], beta_t.ACe[,5], type="o", col="purple")
lines(beta_t.ACe[,1], beta_t.ACe[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y,
c("Intercept", "T(RMRF)", "T(SMB)", "T(HML)", "T(MOM)"),
cex=0.6, col=c("red", "green", "blue", "purple", "brown"),
lty=1, bty="n")

min_y <- min(beta_t.ACv[,2:6])
max_y <- max(beta_t.ACv[,2:6])
plot(beta_t.ACv[,1], beta_t.ACv[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AVW) % Change in T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_t.ACv[,1], beta_t.ACv[,3], type="o", col="green")
lines(beta_t.ACv[,1], beta_t.ACv[,4], type="o", col="blue")
lines(beta_t.ACv[,1], beta_t.ACv[,5], type="o", col="purple")
lines(beta_t.ACv[,1], beta_t.ACv[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y,
c("Intercept", "T(RMRF)", "T(SMB)", "T(HML)", "T(MOM)"),
cex=0.6, col=c("red", "green", "blue", "purple", "brown"),
lty=1, bty="n")

min_y <- min(beta_nwt.ACe[,2:6])
max_y <- max(beta_nwt.ACe[,2:6])
plot(beta_nwt.ACe[,1], beta_nwt.ACe[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AEW) % Change in NW T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,3], type="o", col="green")
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,4], type="o", col="blue")
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,5], type="o", col="purple")
lines(beta_nwt.ACe[,1], beta_nwt.ACe[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y,
c("Intercept", "NWT(RMRF)", "NWT(SMB)", "NWT(HML)", "NWT(MOM)"),
cex=0.6, col=c("red", "green", "blue", "purple", "brown"),
lty=1, bty="n")

min_y <- min(beta_nwt.ACv[,2:6])
max_y <- max(beta_nwt.ACv[,2:6])
plot(beta_nwt.ACv[,1], beta_nwt.ACv[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AVW) % Change in NW T-Stats\nFrom Panel A to Panel C")
axis(1, at=1:10, lab=1:10)
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,3], type="o", col="green")
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,4], type="o", col="blue")
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,5], type="o", col="purple")
lines(beta_nwt.ACv[,1], beta_nwt.ACv[,6], type="o", col="brown")
lines(1:10, rep(0, 10), lty = 3)
legend(10.5, max_y,
c("Intercept", "NWT(RMRF)", "NWT(SMB)", "NWT(HML)", "NWT(MOM)"),
cex=0.6, col=c("red", "green", "blue", "purple", "brown"),
lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### OLS Assumptions 

### Panel B

  ### P-Values of Lilliefors Test for Normality
  normtestpval.B
  
  ### P-values of Non-constant Variance Score Test
  ncvtestpval.B

### Panel C

  ### P-Values of Lilliefors Test for Normality
  normtestpval.C
  
  ### P-values of Non-constant Variance Score Test
  ncvtestpval.C

#############################################################################
### Chapter 4: Portfolio Optimization
#############################################################################

### Mean and SD of Monthly Portfolio Returns

### Plots of Mean and SD

### Fill by column
par(mfcol=c(3,2))
### Expand right side of clipping rect to make room for the legend
par(xpd = T, mar = par()$mar + c(0,0,0,3))

min_y <- min(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
max_y <- max(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
plot(meansd.Ae[,1], meansd.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Mean",
main = "(AEW) Mean of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Be[,1], meansd.Be[,2], type="o", col="red")
lines(meansd.Ce[,1], meansd.Ce[,2], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6,
col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
max_y <- max(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
plot(meansd.Ae[,1], meansd.Ae[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "SD",
main = "(AEW) SD of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Be[,1], meansd.Be[,3], type="o", col="red")
lines(meansd.Ce[,1], meansd.Ce[,3], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6,
col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.ABe[,2:3], meansd.ACe[,2:3])
max_y <- max(meansd.ABe[,2:3], meansd.ACe[,2:3])
plot(meansd.ABe[,1], meansd.ABe[,2], type = "o", col = "skyblue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AEW) % Change from Panel A to Panels B & C")
axis(1, at=1:10, lab=1:10)
lines(meansd.ACe[,1], meansd.ACe[,2], type="o", col="orange")
lines(meansd.ABe[,1], meansd.ABe[,3], type="o", col="purple")
lines(meansd.ACe[,1], meansd.ACe[,3], type="o", col="brown")
legend(10.5, max_y, c("Mean (B)", "Mean (C)", "SD (B)", "SD (C)"), cex=0.6,
col=c("skyblue", "orange", "purple", "brown"), lty=1, bty="n")

min_y <- min(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
max_y <- max(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
plot(meansd.Av[,1], meansd.Av[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "Mean",
main = "(AVW) Mean of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Bv[,1], meansd.Bv[,2], type="o", col="red")
lines(meansd.Cv[,1], meansd.Cv[,2], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6,
col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
max_y <- max(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
plot(meansd.Av[,1], meansd.Av[,3], type = "o", col = "blue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "SD",
main = "(AVW) SD of Monthly Portfolio Return")
axis(1, at=1:10, lab=1:10)
lines(meansd.Bv[,1], meansd.Bv[,3], type="o", col="red")
lines(meansd.Cv[,1], meansd.Cv[,3], type="o", col="green")
legend(10.5, max_y, c("Panel A", "Panel B", "Panel C"), cex=0.6,
col=c("blue", "red", "green"), lty=1, bty="n")

min_y <- min(meansd.ABv[,2:3], meansd.ACv[,2:3])
max_y <- max(meansd.ABv[,2:3], meansd.ACv[,2:3])
plot(meansd.ABv[,1], meansd.ABv[,2], type = "o", col = "skyblue",
ylim = c(min_y,max_y), xlab = "Portfolio", ylab = "%",
main = "(AVW) % Change from Panel A to Panels B & C")
axis(1, at=1:10, lab=1:10)
lines(meansd.ACv[,1], meansd.ACv[,2], type="o", col="orange")
lines(meansd.ABv[,1], meansd.ABv[,3], type="o", col="purple")
lines(meansd.ACv[,1], meansd.ACv[,3], type="o", col="brown")
legend(10.5, max_y, c("Mean (B)", "Mean (C)", "SD (B)", "SD (C)"), cex=0.6,
col=c("skyblue", "orange", "purple", "brown"), lty=1, bty="n")

# Restore default clipping rect
par(mar = c(5, 4, 4, 2) + 0.1)
par(mfrow=c(1,1))

#############################################################################

### Plots of Mean vs SD

par(mfrow=c(3,2))

min_x <- min(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
max_x <- max(meansd.Ae[,3], meansd.Be[,3], meansd.Ce[,3])
min_y <- min(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
max_y <- max(meansd.Ae[,2], meansd.Be[,2], meansd.Ce[,2])
plot(meansd.Ae[,3], meansd.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlim = c(min_x,max_x),
xlab = "SD", ylab = "Mean", main = "(AEW) Mean vs SD")
lines(meansd.Be[,3], meansd.Be[,2], type="o", col="red")
lines(meansd.Ce[,3], meansd.Ce[,2], type="o", col="brown")
legend("topright", c("Panel A", "Panel B", "Panel C"), cex=0.6,
col=c("blue", "red", "brown"), lty=1, bty="n")

min_x <- min(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
max_x <- max(meansd.Av[,3], meansd.Bv[,3], meansd.Cv[,3])
min_y <- min(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
max_y <- max(meansd.Av[,2], meansd.Bv[,2], meansd.Cv[,2])
plot(meansd.Av[,3], meansd.Av[,2], type = "o", col = "green",
ylim = c(min_y,max_y), xlim = c(min_x,max_x),
xlab = "SD", ylab = "Mean", main = "(AVW) Mean vs SD")
lines(meansd.Bv[,3], meansd.Bv[,2], type="o", col="purple")
lines(meansd.Cv[,3], meansd.Cv[,2], type="o", col="orange")
legend("topright", c("Panel A", "Panel B", "Panel C"), cex=0.6,
col=c("green", "purple", "orange"), lty=1, bty="n")

min_x <- min(meansd.Ae[,3], meansd.Av[,3])
max_x <- max(meansd.Ae[,3], meansd.Av[,3])
min_y <- min(meansd.Ae[,2], meansd.Av[,2])
max_y <- max(meansd.Ae[,2], meansd.Av[,2])
plot(meansd.Ae[,3], meansd.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlim = c(min_x,max_x),
xlab = "SD", ylab = "Mean", main = "(Panel A) Mean vs SD")
lines(meansd.Av[,3], meansd.Av[,2], type="o", col="green")
legend("topright", c("AEW", "AVW"), cex=0.6, col=c("blue", "green"),
lty=1, bty="n")

min_x <- min(meansd.Be[,3], meansd.Bv[,3])
max_x <- max(meansd.Be[,3], meansd.Bv[,3])
min_y <- min(meansd.Be[,2], meansd.Bv[,2])
max_y <- max(meansd.Be[,2], meansd.Bv[,2])
plot(meansd.Be[,3], meansd.Be[,2], type = "o", col = "red",
ylim = c(min_y,max_y), xlim = c(min_x,max_x),
xlab = "SD", ylab = "Mean", main = "(Panel B) Mean vs SD")
lines(meansd.Bv[,3], meansd.Bv[,2], type="o", col="purple")
legend("topright", c("AEW", "AVW"), cex=0.6, col=c("red", "purple"),
lty=1, bty="n")

min_x <- min(meansd.Ce[,3], meansd.Cv[,3])
max_x <- max(meansd.Ce[,3], meansd.Cv[,3])
min_y <- min(meansd.Ce[,2], meansd.Cv[,2])
max_y <- max(meansd.Ce[,2], meansd.Cv[,2])
plot(meansd.Ce[,3], meansd.Ce[,2], type = "o", col = "brown",
ylim = c(min_y,max_y), xlim = c(min_x,max_x),
xlab = "SD", ylab = "Mean", main = "(Panel C) Mean vs SD")
lines(meansd.Cv[,3], meansd.Cv[,2], type="o", col="orange")
legend("topright", c("AEW", "AVW"), cex=0.6, col=c("brown", "orange"),
lty=1, bty="n")

min_x <- min(meansd.Ae[,3], meansd.Av[,3],
meansd.Be[,3], meansd.Bv[,3],
meansd.Ce[,3], meansd.Cv[,3])
max_x <- max(meansd.Ae[,3], meansd.Av[,3],
meansd.Be[,3], meansd.Bv[,3],
meansd.Ce[,3], meansd.Cv[,3])
min_y <- min(meansd.Ae[,2], meansd.Av[,2],
meansd.Be[,2], meansd.Bv[,2],
meansd.Ce[,2], meansd.Cv[,2])
max_y <- max(meansd.Ae[,2], meansd.Av[,2],
meansd.Be[,2], meansd.Bv[,2],
meansd.Ce[,2], meansd.Cv[,2])
plot(meansd.Ae[,3], meansd.Ae[,2], type = "o", col = "blue",
ylim = c(min_y,max_y), xlim = c(min_x,max_x),
xlab = "SD", ylab = "Mean", main = "Mean vs SD")
lines(meansd.Av[,3], meansd.Av[,2], type="o", col="green")
lines(meansd.Be[,3], meansd.Be[,2], type = "o", col = "red")
lines(meansd.Bv[,3], meansd.Bv[,2], type="o", col="purple")
lines(meansd.Ce[,3], meansd.Ce[,2], type = "o", col = "brown")
lines(meansd.Cv[,3], meansd.Cv[,2], type="o", col="orange")
legend("topright",
c("AEW (A)", "AVW (A)", "AEW (B)", "AVW (B)", "AEW (C)", "AVW (C)"),
cex=0.6,
col=c("blue", "green", "red", "purple", "brown", "orange"),
lty=1, bty="n")

par(mfrow=c(1,1))

#############################################################################

### Mean-Variance Optimal Portfolio

par(mfrow=c(3,1))

### (AEW) Mean-Variance Optimal Portfolio
plot(effpfT.e[,2], effpfT.e[,1], type = "o", col = "red", 
     xlim = c(min_x.e,max_x.e), xlab = "Standard Deviation",
     ylab = "Expected Monthly Return", 
     main = "(AEW) Mean-Variance Optimal Portfolio")
lines(effpfF.e[,2], effpfF.e[,1], type="o", col="blue")
legend("right", c("Short Allowed", "Long Only"), cex=0.6,
col=c("red", "blue"), lty=1, bty="n")

### (AVW) Mean-Variance Optimal Portfolio
plot(effpfT.v[,2], effpfT.v[,1], type = "o", col = "red",
     xlim = c(min_x.v,max_x.v), xlab = "Standard Deviation",
     ylab = "Expected Monthly Return", 
     main = "(AVW) Mean-Variance Optimal Portfolio")
lines(effpfF.v[,2], effpfF.v[,1], type="o", col="blue")
legend("right", c("Short Allowed", "Long Only"), cex=0.6,
col=c("red", "blue"), lty=1, bty="n")

### (Shorts Allowed) Mean-Variance Optimal Portfolio: AEW vs AVW
plot(effpf.e[,2], effpf.e[,1], type = "o", col = "red",
     xlim = c(min_x.ev,max_x.ev), xlab = "Standard Deviation",
     ylab = "Expected Monthly Return",
     main = "(Shorts Allowed) Mean-Variance Optimal Portfolio: AEW vs AVW")
lines(effpf.v[,2], effpf.v[,1], type="o", col="blue")
legend("right", c("AEW", "AVW"), cex=0.6,
col=c("red", "blue"), lty=1, bty="n")

par(mfrow=c(1,1))

### Mean-Variance Optimal Portfolio: Solution and Correlation

print(
xtable(effpf.e,
caption = "(AEW) Mean-Variance Optimal Portfolio: Solution",
label = "table:t4.1", digits = c(0, rep(3, 2), rep(2, 10)),
display = c('e', rep('f', 12))),
table.placement = "htbp", caption.placement = "top",
include.rownames = FALSE, latex.environment = "center")

print(
xtable(effpf.v,
caption = "(AVW) Mean-Variance Optimal Portfolio: Solution",
label = "table:t4.2", digits = c(0, rep(3, 2), rep(2, 10)),
display = c('e', rep('f', 12))),
table.placement = "htbp", caption.placement = "top",
include.rownames = FALSE, latex.environment = "center")

print(
xtable(cor(effpf.e),
caption = "(AEW) Mean-Variance Optimal Portfolio: Correlation",
label = "table:t4.3", digits = c(0, rep(1, 12)),
display = c('e', rep('f', 12))),
table.placement = "htbp", caption.placement = "top",
latex.environment = "center")

print(
xtable(cor(effpf.v),
caption = "(AVW) Mean-Variance Optimal Portfolio: Correlation",
label = "table:t4.4", digits = c(0, rep(1, 12)),
display = c('e', rep('f', 12))),
table.placement = "htbp", caption.placement = "top",
latex.environment = "center")

### Scatterplot of the Portfolio's Mean, SD, and Weight
pairs(effpfT.e, col="blue")
pairs(effpfT.v, col="blue")

#############################################################################
### Chapter 5: Conclusion
#############################################################################

### None

#############################################################################
### The End of R codes for Chapter
#############################################################################

#############################################################################
### The Beginning of R codes for Appendix
#############################################################################
### Appendix A: Residual Plots
#############################################################################

### (AEW) Portfolio 1
portfolio <- 1
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 1
portfolio <- 1
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 2
portfolio <- 2
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 2
portfolio <- 2
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 3
portfolio <- 3
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 3
portfolio <- 3
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 4
portfolio <- 4
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 4
portfolio <- 4
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 5
portfolio <- 5
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 5
portfolio <- 5
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 6
portfolio <- 6
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 6
portfolio <- 6
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 7
portfolio <- 7
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 7
portfolio <- 7
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 8
portfolio <- 8
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 8
portfolio <- 8
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 9
portfolio <- 9
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 9
portfolio <- 9
residplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 10
portfolio <- 10
residplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 10
portfolio <- 10
residplotFun(data.l, 2, portfolio)

#############################################################################
### Appendix B: Plot of Cook's Distance and Covariance ratio
#############################################################################

### Portfolio 1
toplotFun(data.l, 1, TRUE)

### Portfolio 2
toplotFun(data.l, 2, TRUE)

### Portfolio 3
toplotFun(data.l, 3, TRUE)

### Portfolio 4
toplotFun(data.l, 4, TRUE)

### Portfolio 5
toplotFun(data.l, 5, TRUE)

### Portfolio 6
toplotFun(data.l, 6, TRUE)

### Portfolio 7
toplotFun(data.l, 7, TRUE)

### Portfolio 8
toplotFun(data.l, 8, TRUE)

### Portfolio 9
toplotFun(data.l, 9, TRUE)

### Portfolio 10
toplotFun(data.l, 10, TRUE)

#############################################################################
### Appendix C: Histogram of Cook's Distance and Covariance ratio
#############################################################################

### Portfolio 1
toplotFun(data.l, 1, FALSE)

### Portfolio 2
toplotFun(data.l, 2, FALSE)

### Portfolio 3
toplotFun(data.l, 3, FALSE)

### Portfolio 4
toplotFun(data.l, 4, FALSE)

### Portfolio 5
toplotFun(data.l, 5, FALSE)

### Portfolio 6
toplotFun(data.l, 6, FALSE)

### Portfolio 7
toplotFun(data.l, 7, FALSE)

### Portfolio 8
toplotFun(data.l, 8, FALSE)

### Portfolio 9
toplotFun(data.l, 9, FALSE)

### Portfolio 10
toplotFun(data.l, 10, FALSE)

#############################################################################
### Appendix D: Information for the Excluded Influential Observations
#############################################################################

### (AEW) Portfolio 1
portfolio <- 1
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 1
portfolio <- 1
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 2
portfolio <- 2
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 2
portfolio <- 2
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 3
portfolio <- 3
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 3
portfolio <- 3
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 4
portfolio <- 4
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 4
portfolio <- 4
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 5
portfolio <- 5
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 5
portfolio <- 5
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 6
portfolio <- 6
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 6
portfolio <- 6
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 7
portfolio <- 7
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 7
portfolio <- 7
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 8
portfolio <- 8
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 8
portfolio <- 8
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 9
portfolio <- 9
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 9
portfolio <- 9
exclusioninfoFun(data.v, portfolio)

### (AEW) Portfolio 10
portfolio <- 10
exclusioninfoFun(data.e, portfolio)

### (AVW) Portfolio 10
portfolio <- 10
exclusioninfoFun(data.v, portfolio)

#############################################################################
### Appendix E: Excluded Months By Panel and Weighting Method
#############################################################################

list(
"Panel A's Number of Monthly Observations" = nrow(data.e),
"(AEW) Panel B's Excluded Months" = exclusion.Be,
"(AEW) Panel B's Number of Excluded Months" = length(exclusion.Be),
"(AVW) Panel B's Excluded Months" = exclusion.Bv,
"(AVW) Panel B's Number of Excluded Months" = length(exclusion.Bv),
"(AEW) Panel C's Excluded Months" =
matrix(c(as.character(exclusion.Ce), rep("",3)), ncol = 5, byrow = TRUE),
"(AEW) Panel C's Number of Excluded Months" = length(exclusion.Ce),
"(AVW) Panel C's Excluded Months" =
matrix(c(as.character(exclusion.Cv), rep("",3)), ncol = 5, byrow = TRUE),
"(AVW) Panel C's Number of Excluded Months" = length(exclusion.Cv)
)

#############################################################################
### Appendix F: Plot of MOM vs Monthly Portfolio Return
#############################################################################

### (AEW) Portfolio 1
portfolio <- 1
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 1
portfolio <- 1
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 2
portfolio <- 2
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 2
portfolio <- 2
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 3
portfolio <- 3
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 3
portfolio <- 3
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 4
portfolio <- 4
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 4
portfolio <- 4
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 5
portfolio <- 5
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 5
portfolio <- 5
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 6
portfolio <- 6
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 6
portfolio <- 6
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 7
portfolio <- 7
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 7
portfolio <- 7
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 8
portfolio <- 8
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 8
portfolio <- 8
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 9
portfolio <- 9
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 9
portfolio <- 9
momreturnplotFun(data.l, 2, portfolio)

### (AEW) Portfolio 10
portfolio <- 10
momreturnplotFun(data.l, 1, portfolio)

### (AVW) Portfolio 10
portfolio <- 10
momreturnplotFun(data.l, 2, portfolio)

#############################################################################
### The End of R codes for Appendix
#############################################################################
@

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\printbibliography

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%