In [3]:

```
%pylab inline
import os
from ipypublish import nb_setup
```

Populating the interactive namespace from numpy and matplotlib

In [5]:

```
%load_ext rpy2.ipython
```

In [6]:

```
#%load_ext RWinOut
```

"Walking on water and developing software from a specification are easy if both are frozen" -- Edward V. Berard

In this chapter, we develop some expertise in using the R statistical package. See the manual https://cran.r-project.org/doc/manuals/r-release/R-intro.pdf on the R web site. Work through Appendix A, at least the first page. Also see Grant Farnsworth's document "Econometrics in R": https://cran.r-project.org/doc/contrib/Farnsworth-EconometricsInR.pdf.

There is also a great book that I personally find to be of very high quality, titled "The Art of R Programming" by Norman Matloff.

You can easily install the R programming language, which is a very useful tool for Machine Learning. See: http://en.wikipedia.org/wiki/Machine_learning

Get R from: http://www.r-project.org/ (download and install it).

If you want to use R in IDE mode, download RStudio: http://www.rstudio.com.

Here is quick test to make sure your installation of R is working along with graphics capabilities.

In [7]:

```
%%R
#PLOT HISTOGRAM FROM STANDARD NORMAL RANDOM NUMBERS
x = rnorm(1000000)
hist(x,50)
grid(col="blue",lwd=2)
```

If you want to directly access the system you can issue system commands as follows:

In [8]:

```
%%R
#SYSTEM COMMANDS
#The following command will show the files in the directory which are notebooks.
print(system("ls -lt")) #This command will not work in the notebook.
```

[1] 0

To get started, we need to grab some data. Go to Yahoo! Finance and download some historical data in an Excel spreadsheet, re-sort it into chronological order, then save it as a CSV file. Read the file into R as follows.

In [9]:

```
%%R
#READ IN DATA FROM CSV FILE
data = read.csv("DSTMAA_data/goog.csv",header=TRUE)
print(head(data))
m = length(data)
n = length(data[,1])
print(c("Number of columns = ",m))
print(c("Length of data series = ",n))
```

In [10]:

```
%%R
#REVERSE ORDER THE DATA (Also get some practice with a for loop)
for (j in 1:m) {
data[,j] = rev(data[,j])
}
print(head(data))
stkp = as.matrix(data[,7])
plot(stkp,type="l",col="blue")
grid(lwd=2)
```

We can do the same data set up exercise for financial data using the **quantmod** package.

*Note*: to install a package you can use the drop down menus on Windows and Mac operating systems inside RStudio, and use a package installer on Linux.

You can install R packages from the console using conda: conda install r-

Or issue the following command from the notebook:

(when asked for "selection" enter 60 for California)

In [15]:

```
%%R
#install.packages("quantmod")
```

NULL

Now we move on to using this package for one stock.

In [16]:

```
%%R
#USE THE QUANTMOD PACKAGE TO GET STOCK DATA
library(quantmod)
getSymbols("IBM")
```

/Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Loading required package: xts warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Loading required package: zoo warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Attaching package: ‘zoo’ warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: The following objects are masked from ‘package:base’: as.Date, as.Date.numeric warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Registered S3 method overwritten by 'xts': method from as.zoo.xts zoo warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Loading required package: TTR warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Registered S3 method overwritten by 'quantmod': method from as.zoo.data.frame zoo warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Version 0.4-0 included new data defaults. See ?getSymbols. Learn from a quantmod author: https://www.datacamp.com/courses/importing-and-managing-financial-data-in-r warnings.warn(x, RRuntimeWarning) /Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: ‘getSymbols’ currently uses auto.assign=TRUE by default, but will use auto.assign=FALSE in 0.5-0. You will still be able to use ‘loadSymbols’ to automatically load data. getOption("getSymbols.env") and getOption("getSymbols.auto.assign") will still be checked for alternate defaults. This message is shown once per session and may be disabled by setting options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details. warnings.warn(x, RRuntimeWarning)

[1] "IBM"

In [17]:

```
%%R
chartSeries(IBM)
```

Let's take a quick look at the data.

In [18]:

```
%%R
head(IBM)
```

Extract the dates using pipes (we will see this in more detail later).

In [19]:

```
%%R
library(magrittr)
dts = IBM %>% as.data.frame %>% row.names
dts %>% head %>% print
dts %>% length %>% print
```

[1] "2007-01-03" "2007-01-04" "2007-01-05" "2007-01-08" "2007-01-09" [6] "2007-01-10" [1] 3205

Plot the data.

In [20]:

```
%%R
stkp = as.matrix(IBM$IBM.Adjusted)
rets = diff(log(stkp))
dts = as.Date(dts)
plot(dts,stkp,type="l",col="blue",xlab="Years",ylab="Stock Price of IBM")
grid(lwd=2)
```

Summarize the data.

In [21]:

```
%%R
#DESCRIPTIVE STATS
summary(IBM)
```

Compute risk (volatility).

In [22]:

```
%%R
#STOCK VOLATILITY
sigma_daily = sd(rets)
sigma_annual = sigma_daily*sqrt(252)
print(sigma_annual)
print(c("Sharpe ratio = ",mean(rets)*252/sigma_annual))
```

[1] 0.2217987 [1] "Sharpe ratio = " "0.26146725711648"

We may also use the package to get data for more than one stock.

In [23]:

```
%%R
library(quantmod)
getSymbols(c("GOOG","AAPL","CSCO","IBM"))
```

[1] "GOOG" "AAPL" "CSCO" "IBM"

We now go ahead and concatenate columns of data into one stock data set.

In [24]:

```
%%R
goog = as.numeric(GOOG[,6])
aapl = as.numeric(AAPL[,6])
csco = as.numeric(CSCO[,6])
ibm = as.numeric(IBM[,6])
stkdata = cbind(goog,aapl,csco,ibm)
dim(stkdata)
```

[1] 3205 4

Now, compute daily returns. This time, we do log returns in continuous-time. The mean returns are:

In [25]:

```
%%R
n = dim(stkdata)[1]
rets = log(stkdata[2:n,]/stkdata[1:(n-1),])
colMeans(rets)
```

goog aapl csco ibm 0.0005235354 0.0009525354 0.0002568174 0.0002301313

We can also compute the covariance matrix and correlation matrix:

In [26]:

```
%%R
cv = cov(rets)
print(cv,2)
cr = cor(rets)
print(cr,4)
```

Notice the print command allows you to choose the number of significant digits (in this case 4). Also, as expected the four return time series are positively correlated with each other.

Data frames are the most essential data structure in the R programming language. One may think of a data frame as simply a spreadsheet. In fact you can view it as such with the following command.

In [27]:

```
%%R
#Only works in RStudio
#View(data)
```

NULL

However, data frames in R are much more than mere spreadsheets, which is why Excel will never trump R in the hanlding and analysis of data, except for very small applications on small spreadsheets. One may also think of data frames as databases, and there are many commands that we may use that are database-like, such as joins, merges, filters, selections, etc. Indeed, packages such as **dplyr** and **data.table** are designed to make these operations seamless, and operate efficiently on big data, where the number of observations (rows) are of the order of hundreds of millions.

Data frames can be addressed by column names, so that we do not need to remember column numbers specifically. If you want to find the names of all columns in a data frame, the **names** function does the trick. To address a chosen column, append the column name to the data frame using the "$" connector, as shown below.

In [28]:

```
%%R
#THIS IS A DATA FRAME AND CAN BE REFERENCED BY COLUMN NAMES
print(names(data))
print(head(data$Close))
```

The command printed out the first few observations in the column "Close". All variables and functions in R are "objects", and you are well-served to know the object *type*, because objects have properties and methods apply differently to objects of various types. Therefore, to check an object type, use the **class** function.

In [29]:

```
%%R
class(data)
```

[1] "data.frame"

To obtain descriptive statistics on the data variables in a data frame, the **summary** function is very handy.

In [30]:

```
%%R
#DESCRIPTIVE STATISTICS
summary(data)
```

Let's take a given column of data and perform some transformations on it. We can also plot the data, with some arguments for look and feel, using the **plot** function.

In [31]:

```
%%R
#USING A PARTICULAR COLUMN
stkp = data$Adj.Close
dt = data$Date
print(c("Length of stock series = ",length(stkp)))
#Ln of differenced stk prices gives continuous returns
rets = diff(log(stkp)) #diff() takes first differences
print(c("Length of return series = ",length(rets)))
print(head(rets))
plot(rets,type="l",col="blue")
```

In case you want more descriptive statistics than provided by the **summary** function, then use an appropriate package. We may be interested in the higher-order moments, and we use the **moments** package for this.

In [32]:

```
%%R
print(summary(rets))
```

Compute the daily and annualized standard deviation of returns.

In [33]:

```
%%R
r_sd = sd(rets)
r_sd_annual = r_sd*sqrt(252)
print(c(r_sd,r_sd_annual))
#What if we take the stdev of annualized returns?
print(sd(rets*252))
#Huh?
print(sd(rets*252))/252
print(sd(rets*252))/sqrt(252)
```

[1] 0.01896794 0.30110676 [1] 4.779922 [1] 4.779922 [1] 4.779922 [1] 0.3011068

Notice the interesting use of the **print** function here. The variance is easy as well.

In [34]:

```
%%R
#Variance
r_var = var(rets)
r_var_annual = var(rets)*252
print(c(r_var,r_var_annual))
```

[1] 0.0003597829 0.0906652783

Skewness and kurtosis are key moments that arise in all return distributions. We need a different library in R for these. We use the **moments** library.

Skewness means one tail is fatter than the other (asymmetry). Fatter right (left) tail implies positive (negative) skewness.

\begin{equation} \mbox{Kurtosis} = \frac{E[(X-\mu)^4]}{\sigma^{4}} \end{equation}Kurtosis means both tails are fatter than with a normal distribution.

In [36]:

```
%%R
#HIGHER-ORDER MOMENTS
library(moments)
hist(rets,50)
print(c("Skewness=",skewness(rets)))
print(c("Kurtosis=",kurtosis(rets)))
```

[1] "Skewness=" "-0.57512743445236" [1] "Kurtosis=" "12.6866762557976"

For the normal distribution, skewness is zero, and kurtosis is 3. Kurtosis minus three is denoted "excess kurtosis".

In [38]:

```
%%R
print(skewness(rnorm(1000000)))
print(kurtosis(rnorm(1000000)))
```

[1] 0.00078868 [1] 2.994231

What is the skewness and kurtosis of the stock index (S\&P500)?

Often the original data is in a space delimited file, not a comma separated one, in which case the **read.table** function is appropriate.

In [39]:

```
%%R
#READ IN MORE DATA USING SPACE DELIMITED FILE
data = read.table("DSTMAA_data/markowitzdata.txt",header=TRUE)
print(head(data))
print(c("Length of data series = ",length(data$X.DATE)))
```

We compute covariance and correlation in the data frame.

In [40]:

```
%%R
#COMPUTE COVARIANCE AND CORRELATION
rets = as.data.frame(cbind(data$SUNW,data$MSFT,data$IBM,data$CSCO,data$AMZN))
names(rets) = c("SUNW","MSFT","IBM","CSCO","AMZN")
print(cov(rets))
print(cor(rets))
```

We may redo the example above using a very useful package called **magrittr** which mimics pipes in the Unix operating system. In the code below, we pipe the returns data into the correlation function and then "pipe" the output of that into the print function. This is analogous to issuing the command *print(cor(rets))*.

In [41]:

```
%%R
#Repeat the same process using pipes
library(magrittr)
rets %>% cor %>% print
```

Question: What do you get if you cross a mountain-climber with a mosquito?Answer: Can't be done. You'll be crossing a scaler with a vector.

We will use matrices extensively in modeling, and here we examine the basic commands needed to create and manipulate matrices in R. We create a $4 \times 3$ matrix with random numbers as follows:

In [42]:

```
%%R
x = matrix(rnorm(12),4,3)
print(x)
```

Transposing the matrix, notice that the dimensions are reversed.

In [43]:

```
%%R
print(t(x),3)
```

Of course, it is easy to multiply matrices as long as they conform. By "conform" we mean that when multiplying one matrix by another, the number of columns of the matrix on the left must be equal to the number of rows of the matrix on the right. The resultant matrix that holds the answer of this computation will have the number of rows of the matrix on the left, and the number of columns of the matrix on the right. See the examples below:

In [44]:

```
%%R
print(t(x) %*% x,3)
print(x %*% t(x),3)
```

Here is an example of non-conforming matrices.

In [45]:

```
%%R
#CREATE A RANDOM MATRIX
x = matrix(runif(12),4,3)
print(x)
print(x*2)
print(x+x)
print(t(x) %*% x) #THIS SHOULD BE 3x3
#print(x %*% x) #SHOULD GIVE AN ERROR
```

Taking the inverse of the covariance matrix, we get:

In [46]:

```
%%R
cv_inv = solve(cv)
print(cv_inv,3)
```

Check that the inverse is really so!

In [47]:

```
%%R
print(cv_inv %*% cv,3)
```

It is, the result of multiplying the inverse matrix by the matrix itself results in the identity matrix.

A covariance matrix should be positive definite. Why? What happens if it is not? Checking for this property is easy.

In [49]:

```
%%R
library(corpcor)
is.positive.definite(cv)
```

[1] TRUE

What happens if you compute pairwise covariances from differing lengths of data for each pair?

Let's take the returns data we have and find the inverse.

In [50]:

```
%%R
cv = cov(rets)
print(round(cv,6))
cv_inv = solve(cv) #TAKE THE INVERSE
print(round(cv_inv %*% cv,2)) #CHECK THAT WE GET IDENTITY MATRIX
```

In [51]:

```
%%R
#CHECK IF MATRIX IS POSITIVE DEFINITE (why do we check this?)
library(corpcor)
is.positive.definite(cv)
```

[1] TRUE

Finding roots of nonlinear equations is often required, and R has several packages for this purpose. Here we examine a few examples. Suppose we are given the function
[
(x^2 + y^2 - 1)^3 - x^2 y^3 = 0
]
and for various values of $y$ we wish to solve for the values of $x$. The function we use is called **multiroot** and the use of the function is shown below.

In [53]:

```
%%R
#ROOT SOLVING IN R
library(rootSolve)
fn = function(x,y) {
result = (x^2+y^2-1)^3 - x^2*y^3
}
yy = 1
sol = multiroot(f=fn,start=1,maxiter=10000,rtol=0.000001,atol=0.000001,ctol=0.00001,y=yy)
print(c("solution=",sol$root))
check = fn(sol$root,yy)
print(check)
```

[1] "solution=" "1" [1] 0

Here we demonstrate the use of another function called **uniroot**.

In [54]:

```
%%R
fn = function(x) {
result = 0.065*(x*(1-x))^0.5- 0.05 +0.05*x
}
sol = uniroot.all(f=fn,c(0,1))
print(sol)
check = fn(sol)
print(check)
```

[1] 1.0000000 0.3717627 [1] 0.000000e+00 1.041576e-06

In a *multivariate* linear regression, we have

where $Y \in R^{t \times 1}$, $X \in R^{t \times n}$, and $\beta \in R^{n \times 1}$, and the regression solution is simply equal to $\beta = (X'X)^{-1}(X'Y) \in R^{n \times 1}$.

To get this result we minimize the sum of squared errors.

\begin{eqnarray*} \min_{\beta} e'e &=& (Y - X \cdot \beta)' (Y-X \cdot \beta) \\ &=& Y'(Y-X \cdot \beta) - (X \beta)'\cdot (Y-X \cdot \beta) \\ &=& Y'Y - Y' X \beta - (\beta' X') Y + \beta' X'X \beta \\ &=& Y'Y - Y' X \beta - Y' X \beta + \beta' X'X \beta \\ &=& Y'Y - 2Y' X \beta + \beta' X'X \beta \end{eqnarray*}Note that this expression is a scalar.

Differentiating w.r.t. $\beta'$ gives the following f.o.c:

\begin{eqnarray*} - 2 X'Y + 2 X'X \beta&=& {\bf 0} \\ & \Longrightarrow & \\ \beta &=& (X'X)^{-1} (X'Y) \end{eqnarray*}There is another useful expression for each individual $\beta_i = \frac{Cov(X_i,Y)}{Var(X_i)}$. You should compute this and check that each coefficient in the regression is indeed equal to the $\beta_i$ from this calculation.

*Example*: We run a stock return regression to exemplify the algebra above.

In [55]:

```
%%R
data = read.table("DSTMAA_data/markowitzdata.txt",header=TRUE) #THESE DATA ARE RETURNS
print(names(data)) #THIS IS A DATA FRAME (important construct in R)
head(data)
```

In [56]:

```
%%R
#RUN A MULTIVARIATE REGRESSION ON STOCK DATA
Y = as.matrix(data$SUNW)
X = as.matrix(data[,3:6])
res = lm(Y~X)
summary(res)
```

Now we can cross-check the regression using the algebraic solution for the regression coefficients.

In [57]:

```
%%R
#CHECK THE REGRESSION
n = length(Y)
X = cbind(matrix(1,n,1),X)
b = solve(t(X) %*% X) %*% (t(X) %*% Y)
print(b)
```

[,1] -0.0007256342 MSFT 0.1382312148 IBM 0.3791500328 CSCO 0.5769097262 AMZN 0.0324898716

*Example*: As a second example, we take data on basketball teams in a cross-section, and try to explain their performance using team statistics. Here is a simple regression run on some data from the 2005-06 NCAA basketball season for the March madness stats. The data is stored in a space-delimited file called **ncaa.txt**. We use the metric of performance to be the number of games played, with more successful teams playing more playoff games, and then try to see what variables explain it best. We apply a simple linear regression that uses the R command **lm**, which stands for "linear model".

In [58]:

```
%%R
#REGRESSION ON NCAA BASKETBALL PLAYOFF DATA
ncaa = read.table("DSTMAA_data/ncaa.txt",header=TRUE)
print(head(ncaa))
y = ncaa[3]
y = as.matrix(y)
x = ncaa[4:14]
x = as.matrix(x)
```

In [59]:

```
%%R
fm = lm(y~x)
res = summary(fm)
res
```

An alternative specification of regression using data frames is somewhat easier to implement.

In [60]:

```
%%R
#CREATING DATA FRAMES
ncaa_data_frame = data.frame(y=as.matrix(ncaa[3]),x=as.matrix(ncaa[4:14]))
fm = lm(y~x,data=ncaa_data_frame)
summary(fm)
```

In a regression, we estimate the **coefficients** on each explanatory variable (these are the Estimates in the regression above). In addition, we also estimate the standard deviation of the coefficient value, which implies the range around the mean value. This is called the standard error of the coefficient. We are interested in making sure that the coefficient value $b$ is not zero in the statistical sense, usually taken to mean that it is at least 2 standard deviations away from zero. That is, we want the coefficient value $b$ divided by its standard deviation $\sigma_b$ to be at least 2. This is called the **t-statistic** or t-value, shown in the regression above.

The t-statistic is the number of standard deviations the coefficient is away from zero. It implies a **p-value**, which is a probability that the coefficient is equal to zero. So, we want the p-values to be small. We see in the above regression [Pr(>|t|)] that the coefficients that are statistically significant have small p-values and large absolute values of t-statistics. It is intuitive that when the t-statistic is large (negative or positive) it means that the coefficient is far away from zero and using the standard normal distribution we can calculate the probability left in the tails. So if the t-statistic is (say) 2.843, it means that there is only 0.006375 probability remaining in the right tail to the right of the t-statistic value.

For a more detailed discussion see this excellent article in the Scientific American (2019); pdf.

The linear regression is fit by minimizing the sum of squared errors, but the same concept may also be applied to a nonlinear regression as well. So we might have:

$$ y_i = f(x_{i1},x_{i2},...,x_{ip}) + \epsilon_i, \quad i=1,2,...,n $$which describes a data set that has $n$ rows and $p$ columns, which are the standard variables for the number of rows and columns. Note that the error term (residual) is $\epsilon_i$.

The regression will have $(p+1)$ coefficients, i.e., ${\bf b} = \{b_0,b_1,b_2,...,b_p\}$, and ${\bf x}_i = \{x_{i1},x_{i2},...,x_{ip}\}$. The model is fit by minimizing the sum of squared residuals, i.e.,

$$ \min_{\bf b} \sum_{i=1}^n \epsilon_i^2 $$We define the following:

- Sum of squared residuals (errors): $SSE = \sum_{i=1}^n \epsilon_i^2$, with degrees of freedom $DFE= n-p$.
- Total sum of squares: $SST = \sum_{i=1}^n (y_i - {\bar y})^2$, where ${\bar y}$ is the mean of $y$. Degrees of freedom are $DFT = n-1$.
- Regression (model) sum of squares: $SSM = \sum_{i=1}^n (f({\bf x}_i) - {\bar y})^2$; with degrees of freedom $DFM = p-1$.
- Note that $SST = SSM + SSE$.
- Check that $DFT = DFM + DFE$.

The $R$-squared of the regression is

$$ R^2 = \left( 1 - \frac{SSE}{SST} \right) \quad \in (0,1) $$The $F$-statistic in the regression is what tells us if the RHS variables comprise a model that explains the LHS variable sufficiently. Do the RHS variables offer more of an explanation that simply assuming that the mean value of $y$ is the best prediction? The null hypothesis we care about is

- $H_0$: $b_k = 0, k=0,1,2,...,p$, versus an alternate hypothesis of
- $H_1$: $b_k \neq 0$ for at least one $k$.

To test this the $F$-statistic is computed as the following ratio:

$$ F = \frac{\mbox{Explained variance}}{\mbox{Unexplained variance}} = \frac{SSM/DFM}{SSE/DFE} = \frac{MSM}{MSE} $$where $MSM$ is the mean squared model error, and $MSE$ is mean squared error.

Now let's relate this to $R^2$. First, we find an approximation for the $R^2$.

$$ R^2 = 1 - \frac{SSE}{SST} \\ = 1 - \frac{SSE/n}{SST/n} \\ \approx 1 - \frac{MSE}{MST} \\ = \frac{MST-MSE}{MST} \\ = \frac{MSM}{MST} $$The $R^2$ of a regression that has no RHS variables is zero, and of course $MSM=0$. In such a regression $MST = MSE$. So the expression above becomes:

$$ R^2_{p=0} = \frac{MSM}{MST} = 0 $$We can also see with some manipulation, that $R^2$ is related to $F$ (approximately, assuming large $n$).

$$ R^2 + \frac{1}{F+1}=1 \quad \mbox{or} \quad 1+F = \frac{1}{1-R^2} $$Check to see that when $R^2=0$, then $F=0$.

We can further check the formulae with a numerical example, by creating some sample data.

In [61]:

```
%%R
x = matrix(runif(300),100,3)
y = 5 + 4*x[,1] + 3*x[,2] + 2*x[,3] + rnorm(100)
y = as.matrix(y)
res = lm(y~x)
print(summary(res))
```

In [62]:

```
%%R
e = res$residuals
SSE = sum(e^2)
SST = sum((y-mean(y))^2)
SSM = SST - SSE
print(c(SSE,SSM,SST))
R2 = 1 - SSE/SST
print(R2)
n = dim(x)[1]
p = dim(x)[2]+1
MSE = SSE/(n-p)
MSM = SSM/(p-1)
MST = SST/(n-1)
print(c(n,p,MSE,MSM,MST))
Fstat = MSM/MSE
print(Fstat)
```

We can also compare two regressions, say one with 5 RHS variables with one that has only 3 of those five to see whether the additional two variables has any extra value. The ratio of the two $MSM$ values from the first and second regressions is also a $F$-statistic that may be tested for it to be large enough.

Note that if the residuals $\epsilon$ are assumed to be normally distributed, then squared residuals are distributed as per the chi-square ($\chi^2$) distribution. Further, the sum of residuals is distributed normal and the sum of squared residuals is distributed $\chi^2$. And finally, the ratio of two $\chi^2$ variables is $F$-distributed, which is why we call it the $F$-statistic, it is the ratio of two sums of squared errors.

Underlying the analyses of the regression model above is an assumption that the error term $\epsilon$ is independent of the $x$ variables. This assumption ensures that the regression coefficient $\beta$ is unbiased. To see this in the simplest way, consider the univariate regression

$$ y = \beta x + \epsilon $$We have seen earlier that the coefficient $\beta$ is given by

$$ \frac{Cov(x,y)}{Var(x)} = \frac{Cov(x,\beta x + \epsilon)}{Cov(x,x)} = \beta + \frac{Cov(x,\epsilon)}{Cov(x,x)} $$This little piece of statistical math shows that this is biased if there is correlation between $x$ and $\epsilon$.

One way in which the coefficient is biased is if there is a missing variable in the regression that has an effect on both $x$ and $y$, which then injects correlation between $x$ and $\epsilon$. If there is a missing variable that impacts $y$ and not $x$, then it is just fine, after all, every regression has missing variables, else there would be no residual (error) term. Hopefully, there is some idea of how the missing variable impacts both $x$ and $y$ (direction, and if possible sign). Then at least one might have a sense of the direction of bias in the regression coefficient.

Simple linear regression assumes that the standard error of the residuals is the same for all observations. Many regressions suffer from the failure of this condition. The word for this is "heteroskedastic" errors. "Hetero" means different, and "skedastic" means dependent on type.

We can first test for the presence of heteroskedasticity using a standard Breusch-Pagan test available in R. This resides in the **lmtest** package which is loaded in before running the test.

In [64]:

```
%%R
ncaa = read.table("DSTMAA_data/ncaa.txt",header=TRUE)
y = as.matrix(ncaa[3])
x = as.matrix(ncaa[4:14])
result = lm(y~x)
library(lmtest)
bptest(result)
```

studentized Breusch-Pagan test data: result BP = 15.538, df = 11, p-value = 0.1592

We can see that there is very little evidence of heteroskedasticity in the standard errors as the $p$-value is not small. However, lets go ahead and correct the t-statistics for heteroskedasticity as follows, using the **hccm** function. The **hccm** stands for heteroskedasticity corrected covariance matrix.

In [65]:

```
%%R
wuns = matrix(1,64,1)
z = cbind(wuns,x)
b = solve(t(z) %*% z) %*% (t(z) %*% y)
result = lm(y~x)
library(car)
vb = hccm(result)
stdb = sqrt(diag(vb))
tstats = b/stdb
print(tstats)
```

/Users/srdas/anaconda3/lib/python3.7/site-packages/rpy2/rinterface/__init__.py:146: RRuntimeWarning: Loading required package: carData warnings.warn(x, RRuntimeWarning)

GMS -2.68006069 PTS -0.38212818 REB 2.38342637 AST -0.40848721 TO -0.28709450 A.T 0.65632053 STL 2.13627108 BLK 0.09548606 PF -0.68036944 FG 3.52193532 FT 2.35677255 X3P 1.23897636

Compare these to the t-statistics in the original model

In [66]:

```
%%R
summary(result)
```

It is apparent that when corrected for heteroskedasticity, the t-statistics in the regression are lower, and also render some of the previously significant coefficients insignificant.

When data is autocorrelated, i.e., has dependence in time, not accounting for this issue results in unnecessarily high statistical significance (in terms of inflated t-statistics). Intuitively, this is because observations are treated as independent when actually they are correlated in time, and therefore, the true number of observations is effectively less.

Consider a finance application. In efficient markets, the correlation of stock returns from one period to the next should be close to zero. We use the returns on Google stock as an example. First, read in the data.

In [67]:

```
%%R
data = read.csv("DSTMAA_data/goog.csv",header=TRUE)
head(data)
```

Next, create the returns time series.

In [68]:

```
%%R
n = length(data$Close)
stkp = rev(data$Adj.Close)
rets = as.matrix(log(stkp[2:n]/stkp[1:(n-1)]))
n = length(rets)
plot(rets,type="l",col="blue")
print(n)
```

[1] 3606