Open Source: R Programming and Statistical Analysis {#IntroductoryRprogamming}

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

Got R?

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.

System Commands

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

Loading Data

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.

Getting External Stock Data

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)

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

Let's take a quick look at the data.

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

Plot the data.

Summarize the data.

Compute risk (volatility).

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

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

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

We can also compute the covariance matrix and correlation matrix:

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

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.

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.

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.

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

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 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.

Compute the daily and annualized standard deviation of returns.

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

Higher-Order Moments

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.

\begin{equation} \mbox{Skewness} = \frac{E[(X-\mu)^3]}{\sigma^{3}} \end{equation}

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.

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

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

Reading space delimited files

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

We compute covariance and correlation in the data frame.

Pipes with magrittr

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)).

Matrices

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:

Transposing the matrix, notice that the dimensions are reversed.

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:

Here is an example of non-conforming matrices.

Taking the inverse of the covariance matrix, we get:

Check that the inverse is really so!

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.

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.

Root Finding

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.

Here we demonstrate the use of another function called uniroot.

Regression

In a multivariate linear regression, we have

\begin{equation} Y = X \cdot \beta + e \end{equation}

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.

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

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".

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

P-Values, t-statistics

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.

Parts of a regression

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:

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

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.

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.

Bias in regression coefficients

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.

Heteroskedasticity

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.

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.

Compare these to the t-statistics in the original model

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.

Auto-Regressive Models

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.

Next, create the returns time series.