Final project

Contents

1 Basic Commands 2

2 Working with a Data Set 3

3 Basic Linear Regression 4

4 Ridge Regression 4

5 LASSO 5

6 Cross Validation: Ridge vs Lasso 5

7 Forward and backward step wise selection 5

8 Cross Validation: Ridge vs LASSO vs Stepwise 6

9 Polynomials 7

10 Splines 8

11 Natural Cubic Spline 9

12 Generalized Additive Models 9

13 Saving Predictions 9

1

1 Basic Commands

# List of basic commands from the discussion section. ## Clears global enviroment rm(list = ls(all.names = TRUE))

## Define a vector my_first_vector <- c(1,2,3)

## Define a vector with repeated elements my_second_vector <- rep(1, 10)

## Combine the two vectors my_third_vector <- c(my_second_vector, my_first_vector)

## Square the vector print(my_third_vector**2)

## [1] 1 1 1 1 1 1 1 1 1 1 1 4 9

## Define a matrix A <- matrix(c(1:9), nrow = 3, ncol = 3) B <- matrix(c(1:6), nrow = 3, ncol = 3 )

## Transpose of a matrix B_t <- t(B)

## Matrix multiplication A%*%B #Note that the we need to use %*%

## [,1] [,2] [,3] ## [1,] 30 66 30 ## [2,] 36 81 36 ## [3,] 42 96 42

A*B #Component wise multiplication: yields different results

## [,1] [,2] [,3] ## [1,] 1 16 7 ## [2,] 4 25 16 ## [3,] 9 36 27

## Fill a vector with missing values C <- rep(NA, 10) #Useful for loops. Here you will store the results.

## Fill a matrix with missing values D <- matrix(NA, nrow = 2, ncol = 2)

## Add columns to a matrix cbind(A, c(10,11,12))

## [,1] [,2] [,3] [,4] ## [1,] 1 4 7 10 ## [2,] 2 5 8 11 ## [3,] 3 6 9 12

2

## Add rows to a matrix rbind(A, c(10,11,12))

## [,1] [,2] [,3] ## [1,] 1 4 7 ## [2,] 2 5 8 ## [3,] 3 6 9 ## [4,] 10 11 12

2 Working with a Data Set

## Working with a data set: once you have downloaded your data set to your computer, #set the working directory to be the folder where you saved the downloaded data set.

## Set working directoty setwd("/Users/julianmartineziriarte/OneDrive - UC San Diego/Econ 178/Final Project") #Set #the working directory. Here R will look for your file

## Open a data set: data_tr <- read.table("data.txt", header = TRUE, sep = "\t", dec = ".")[,-1] #Load the data #set into R

## Add columns to a data set: data_tr2 <- cbind(data_tr, data_tr$age**2) #Here we added age^2 to our data set.

## Inspect the data set dim(data_tr) #dimension: number of observations and variables.

## [1] 9915 44

names(data_tr) #names of variables: this is useful to call the variables

## [1] "ira" "a401" "hval" "hmort" "hequity" "nifa" "net_nifa" ## [8] "tfa" "net_tfa" "tfa_he" "tw" "age" "inc" "fsize" ## [15] "educ" "db" "marr" "male" "twoearn" "dum91" "e401" ## [22] "p401" "pira" "nohs" "hs" "smcol" "col" "icat" ## [29] "ecat" "zhat" "net_n401" "hown" "i1" "i2" "i3" ## [36] "i4" "i5" "i6" "i7" "a1" "a2" "a3" ## [43] "a4" "a5"

summary (data_tr$tw) #note that to work with a particular variable we #call it data_tr$ + name of variable.

## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -502302 3292 25100 63817 81488 2029910

pairs(data_tr[,c(1,2,3)]) #dependence between the first three covariates

hist(data_tr$tw, breaks=100, main="histogram of total wealth", xlab="total wealth") #Histogram

plot(data_tr$age, data_tr$tw, main="plot of age vs total wealth", xlab="age", ylab="total wealth") #Scatterplot

#Dealing with dummies: use the model.matrix command

3

data(Wage, package="ISLR") new_X <- model.matrix(~. -1, data = Wage[, c(1:9)]) # We transformed the matrix #into a set of 0/1 dummy variables. Now every component is numeric!

3 Basic Linear Regression

reg <- lm(tw ~ age, data=data_tr) #Note we need to specify where the data comes from.

plot(reg) #These multiple plots will allow you to check for outliers #and high leverage observations: see section 3.3.3 of the book.

plot(tw ~ age, data=data_tr, col="red", pch=19, cex=1, main="Age VS Total Wealth") abline(reg, col="blue", lwd=2) #Plots the fitted line

reg2 <- lm(tw ~ age+hmort, data=data_tr) ##Multiple linear regression: you sum the variables reg3 <- lm(tw ~ 1, data=data_tr) #Regression on just a constant reg4 <- lm(tw ~ ., data=data_tr) #Regression on the whole data set: #perfect collinear variables are dropped automatically.

ols_pred <- predict(reg2, newdata = data_tr) #VERY IMPORTANT: computes #predictions of your model on the training data set. You could predict on a different data set. Here we use the same one MSE <- mean((data_tr$tw - ols_pred)**2) # Computes mean square error of your predictions.

4 Ridge Regression

#Install and load package. Please read its description if you have #any doubt: https://cran.r-project.org/web/packages/glmnet/glmnet.pdf #install.packages('glmnet') library(glmnet)

#Choose sequence of lambdas lambdas.rr <- exp(seq(-1, 10, length = 50)) #You have to specify your own sequence.

#Convert your data to matrix: this is very important. Without this #transformation your code will not run. y <- data_tr$tw X <- as.matrix(data_tr[,-1])

ridge_cv <- cv.glmnet(x = X, y = y, lambda = lambdas.rr,alpha = 0) #alpha=0 for Ridge penalty ridge_cv$lambda.min #value of lambda that gives minimum mean cross-validated error

## [1] 0.3678794

#IMPORTANT: Here is how you predict using the cross-valited lambda: ridge_pred <- predict(ridge_cv, newx=X, s = ridge_cv$lambda.min ) #Predict on the

4

#training data set using lambdamin.

5 LASSO

#Choose sequence of lambdas: need not the same as in ridge. lambdas.rr <- exp(seq(-1, 10, length = 50)) #You have to specify your own sequence.

lasso_cv <- cv.glmnet(x = X, y = y, lambda = lambdas.rr,alpha = 1) #alpha=1 for LASSO penalty lasso_cv$lambda.min #value of lambda that gives minimum mean cross-validated error

## [1] 0.3678794

#IMPORTANT: Here is how you predict using the cross-validated lambda: #Predict on the training data set using lambdamin. lasso_pred <- predict(lasso_cv, newx=X, s = lasso_cv$lambda.min ) #Predict on the #training data set using lambdamin.

6 Cross Validation: Ridge vs Lasso

## Fast way to do CV to choose between Lasso and Ridge. rm(list = ls(all.names = TRUE))

data(Wage, package="ISLR") new_X <- model.matrix(~. -1, data = Wage[, c(1:9)])

### Run lasso and ridge regression library(glmnet) my_lasso <- cv.glmnet(y = Wage$logwage, x = as.matrix(new_X), alpha = 1) my_ridge <- cv.glmnet(y = Wage$logwage, x = as.matrix(new_X), alpha = 0)

### Short-cuts: looking at MSPE

## Return the optimal MSPE: we choose the minimum MSPE. c(min(my_lasso$cvm), min(my_ridge$cvm))

## [1] 0.07811426 0.07858603

7 Forward and backward step wise selection

#install.packages('MASS') library(MASS)

# Define the full and null model

5

full <- lm(tw ~ ., data=data_tr) null <- lm(tw ~ 1, data=data_tr)

# forward stepwise - AIC a <- stepAIC(null, scope=list(lower=null, upper=full), direction='forward')

# backward stepwise - AIC b <- stepAIC(full, scope=list(lower=null, upper=full), direction='backward')

# Look at the coefficients to know which variables were selected coef(a)

## (Intercept) tfa_he inc ira twoearn hval ## -8.141847e+03 1.004772e+00 3.348468e-01 3.364930e-01 -5.873304e+03 6.982005e-02 ## nifa net_nifa age hown pira icat ## 4.478408e-01 -4.121149e-01 2.493898e+02 -4.212455e+03 5.372433e+03 -2.087495e+03 ## a3 marr male p401 ## -2.662806e+03 3.911915e+03 3.030529e+03 -2.512205e+03

coef(b)

## (Intercept) ira hval hmort nifa net_nifa ## 6894.1531572 0.3260869 1.0737695 -1.0022450 -0.5692419 0.5928264 ## tfa inc marr male twoearn p401 ## 1.0167355 0.3317287 3915.2788092 3009.1442532 -5890.9702352 -2656.0358267 ## pira icat hown a1 a2 a3 ## 5350.6283321 -2055.4088948 -4181.7142808 -7879.8155593 -7460.7839632 -7806.3235251 ## a4 ## -2871.5036874

# And finally we predict (you could predict on a different data set, #here we use the same one) pr.stepwise_backward <- predict(b, newdata=data_tr) pr.stepwise_forward <- predict(a, newdata=data_tr)

8 Cross Validation: Ridge vs LASSO vs Stepwise

# We do cross validation to choose among ridge, lasso, stepwise backward and stepwise forward.

n <- length(y) #number of observations k <- 10 #number of folds ii <- sample(rep(1:k, length= n)) #randomly assign the observations to each of the 10 folds pr.stepwise_backward <- pr.stepwise_forward <- pr.lasso <- pr.ridge <- rep(NA, length(y)) # prepare the vectors on which to store the values of the cross-validated predictions

#If you have doubts about the loop, the best way to understand #what a loop is doing is: create "j<- 1" and run the loop manually line by line. for (j in 1:k){

hold <- (ii == j)

6

train <- (ii != j) ## Stepwise full <- lm(tw ~ ., data=data_tr[train,]) null <- lm(tw ~ 1, data=data_tr[train,]) a <- stepAIC(null, scope=list(lower=null, upper=full), direction='forward') # backward stepwise - AIC b <- stepAIC(full, scope=list(lower=null, upper=full), direction='backward') pr.stepwise_backward[hold] <- predict(b, newdata=data_tr[hold,]) pr.stepwise_forward[hold] <- predict(a, newdata=data_tr[hold,])

## Do with lass (we use X and y defined above) xx.tr <- X[train,] y.tr <- y[train] xx.te <- X[hold,] ridge.cv <- cv.glmnet(x=as.matrix(xx.tr), y=y.tr, nfolds=k, alpha=0) #glmnet chooses #its own sequence if no lambda sequence is provided. lasso.cv <- cv.glmnet(x=as.matrix(xx.tr), y=y.tr, nfolds=k, alpha=1) #glmnet chooses #its own sequence if no lambda sequence is provided. pr.lasso[hold] <- predict(lasso.cv, newx=xx.te) pr.ridge[hold] <- predict(ridge.cv, newx=xx.te)

}

mspe_step_backward <- mean((pr.stepwise_backward-y)^2) mspe_step_forward <- mean((pr.stepwise_forward-y)^2) mspe.Lasso <- mean((pr.lasso-y)^2) mspe.ridge <- mean((pr.ridge-y)^2)

9 Polynomials

data(lidar, package="SemiPar")

# We use the command "poly" to generate polynomials of a desired order. pol_1 <- lm(logratio ~ poly(range,1), data=lidar) #Does a poynomial regression #with a polynomial of order 1 pol_3 <- lm(logratio ~ poly(range,3), data=lidar) #Does a poynomial regression #with a polynomial of order 3 pol_5 <- lm(logratio ~ poly(range,5), data=lidar) #Does a poynomial regression #with a polynomial of order 5 pol_20 <- lm(logratio ~ poly(range,20), data=lidar) #Does a poynomial regression #with a polynomial of order 20

# We can plot all the polynomials in the same graph to see how they differ plot(cbind(lidar$range, lidar$logratio), xlab="x", ylab="y", main="Lidar", pch=19, col="red") lines(cbind(lidar$range,predict(pol_1)), col="green", lwd=2) lines(cbind(lidar$range,predict(pol_3)), col="blue", lwd=2) lines(cbind(lidar$range,predict(pol_5)), col="gray", lwd=2) lines(cbind(lidar$range,predict(pol_20)), col="black", lwd=2) #Notice the slight #overfit that appears here. Be careful when using high degree polynomials as numerical issues appear.

# Cross Validation to select the degree of the polynomial.

7

y <- lidar$logratio n <- length(y) k <- 10 ii <- sample(rep(1:k, length= n)) p <- 20 #max number of polynomials

pol_cv <- matrix(data=NA,nrow=length(y),ncol=p) mse_cv <- rep(NA, p)

for (s in 1:p){

for (j in 1:k){

hold <- (ii == j) train <- (ii != j) pol <- lm(logratio ~ poly(range,s), data=lidar[train,]) pol_cv[hold,s] <- predict(pol, newdata=lidar[hold,])

} mse_cv[s] <- mean((y - pol_cv[,s])**2)

}

plot(mse_cv)

10 Splines

# To understand the formula for the piecewise polynomials, please see #the explanation in: https://github.com/dviviano/ECON178_TA/blob/master/TA_lectures/Lecture7.ipynb

## Linear spline

kn<- as.numeric( quantile(lidar$range, (1:4)/5) )## Prepare the knots

# Prepare the matrix of covariates x <- matrix(0, dim(lidar)[1], length(kn)+1) for(j in 1:length(kn)) {

x[,j] <- pmax(lidar$range-kn[j], 0) } x[, length(kn)+1] <- lidar$range

# Fit the regression model ppm <- lm(lidar$logratio ~ x) plot(logratio~range, data=lidar, pch=19, col='red', cex=1, main="lidar", xlab="x", ylab="y") lines(predict(ppm)[order(range)]~sort(range), data=lidar, lwd=2, col='blue')

## Cubic spline kn <- as.numeric( quantile(lidar$range, (1:4)/5) )

# prepare the matrix of covariates x <- matrix(0, dim(lidar)[1], length(kn)+3) for(j in 1:length(kn)) {

x[,j] <- pmax(lidar$range-kn[j], 0)^3 ## <<<---- Note: here you put the degree

8

} x[, length(kn)+1] <- lidar$range x[, length(kn)+2] <- lidar$range^2 x[, length(kn)+3] <- lidar$range^3 ## <<------ Note: also here (see formula)

# Fit the regression model ppm <- lm(lidar$logratio ~ x) plot(logratio~range, data=lidar, pch=19, col='red', cex=1, main="lidar", xlab="x", ylab="y") lines(predict(ppm)[order(range)]~sort(range), data=lidar, lwd=2, col='blue')

11 Natural Cubic Spline

library(splines) ## Construct the splines (note ns() create a basis expansion) reg3.ns <- lm(logratio ~ ns(x=range, df=4), data=lidar) #Four degrees of freedom equals to #three interior knots. reg4.ns <- lm(logratio ~ ns(x=range, df=10), data=lidar) #Ten degrees of freedom

plot(cbind(lidar$range, lidar$logratio), xlab="x", ylab="y", main="Lidar", pch=19, col="red") lines(cbind(lidar$range,predict(reg3.ns )), col="green", lwd=2) #natural cubic spline lines(cbind(lidar$range,predict(pol_3)), col="blue", lwd=2) #third degree polynomial: #they look similar, but they are not the same.

## The degrees of freedom must be chosen with cross validation.

12 Generalized Additive Models

#install.packages("ISLR") library(ISLR) data(Wage, package="ISLR")

#For simple generalized dditive models, you can just a linear regression adding the #transformation of each variable. gam1 <- lm(wage ~ ns(year ,4)+ns(age ,5)+education ,data=Wage) #This fits a natural #cubic spline for each variable: year with 4 df and and age with 5 df.

13 Saving Predictions

We have uploaded a file called

"data_for_prediction.txt"

to the "final project" folder on Canvas. This file contains 1982 observations for every variable in the data set

"data_tr.txt"

except for the Total Wealth ("tw") variable.

9

This is what you have to do: 1) Using the data in

"data_tr.txt"

you select a model for prediction. Say for example you settled on a Lasso model (this is a very simplistic example, we hope you do better in your final project):

my_model <- cv.glmnet(y = data_tr$tw, x = as.matrix(data_tr[, c(2:18)], alpha = 1))

2) Using the data in

"data_for_prediction.txt"

and your model, you are going to predict the variable Total Wealth. For example, your code would look something like this

my_predictions <- predict(my_model, newx = as.matrix(data_te[,c(2:18)]))

This generates a vector of 1982 predictions for Total Wealth. You must save this using the following command (don’t forget to set a working directory)

write.table(my_predictions, file = 'my_predictions.txt')

3) Now, we are going to compare your 1982 predictions for Total Wealth stored in

"my_predictions.txt"

to the real value of Total Wealth that we have kept from you. The idea is that if your model works well, it will correctly in this test data set. The better your model predict, the more bonus points you will get.

It is very important that you submit the file

"my_predictions.txt"

that contains your 1982 predictions in order to be eligible for bonus points.

10