辅导 program、R程序设计讲解
Problem Set 1
XXX
Due: 10/10/2024
Introduction
This document was produced by R using RMarkdown. To complete this weeks assignment, we will ask you to
complete a series of analytical and coding exercises. The Analytical Exercises require no coding, whereas
the Coding Exercises require you to use R. The nice thing about RMarkdown is you can do both your
analytical and coding exercises in the same document. For each part in the Coding Exercises, we provide
an empty space of code chunks (area highlighted in grey with a header of the form “## ~~ Problem XX ~~
##”“)
To ease your introduction into R, Problem 2 is a short tutorial into the R programming environment.
Hopefully, you have already downloaded RStudio on your computer. If not, please go do that now. You can
download the latest version at (https://www.rstudio.com/products/rstudio/download/).
Once you have downloaded RStudio, you will be able to open the R markdown script (assignment_1.rmd
ffle) that created this assignment_1.pdf ffle. We ask you to ffll in your code in the code chunk sections (the
areas highlighted in gray bounded by “‘ marks) in the .rmd ffle in each of the subparts of the questions. You
will see that you can include LaTex code in these document. You are not required to use LaTex to do the
analytical exercises (i.e. those without coding), but it is good practice to improve your LaTex skills.
In order to compile a markdown document (i.e. turn your code into a pdf ffle), you must have a version of LaTex
downloaded on your computer. I suggest you download MikTex (https://miktex.org/howto/install-miktex).
If at any time you are confused about R, or not sure what a command does or additional arguments available
for each command, there are two tried and true methods to help resolve this issue. In the R console, you can
use the help command, where you supply the name of the command you are confused about. Alternatively,
google is your friend.
Installing MikTex R Markdown was recently updated, and this update has issues with missing packages
with MikTex. In order to deal with these issues, you need to allow MikTex to install any missing packages
without asking you ffrst. To do so, when you are installing MikTex, in the ‘Settings’ screen, it asks Install
missing packages on-the-ffy. Please select Yes in this screen. If you have already installed MikTex, you
can go to the MikTex console -> Settings and the same box appears in that screen.
Packages to Install
Each week, we will list the packages that you need to install into R in order for you to complete the assignments.
This also allows you to know a nice resource to view which packages you have learned throughout this course.
The packages used this week are
• stats
• ggplot2 (optional)
1Code Setup
## This is a code chunk: it is outlined in grey and has R code inside of it
## Note: you can change what is shown in the final .pdf document using arguments
## inside the curly braces at the top {r, comment='\t\t'}. For example, you
## can turn off print statements showing in the .pdf by adding 'echo=False'
## i.e. changing the header to {r, comment='\t\t',echo=FALSE}
## ~~~~~~~~~~~~~~ CODE SETUP ~~~~~~~~~~~~~~ ##
# ~ This bit of code will be hidden after Problem Set 1 ~
#
# This section sets up the correct directory structure so that
# the working directory for your code is always in the data folder
# Retrieve the code working directory
#script_dir = dirname(sys.frame(1)$ofile)
initial_options <- commandArgs(trailingOnly = FALSE)
render_command <- initial_options[grep('render',initial_options)]
script_name <- gsub("'", "",
regmatches(render_command,
gregexpr("'([ˆ']*)'",
render_command))[[1]][1])
# Determine OS (backslash versus forward slash directory system)
sep_vals = c(length(grep('\\\\',script_name))>0,length(grep('/',script_name))>0)
file_sep = c("\\\\","/")[sep_vals]
# Get data directory
split_str = strsplit(script_name,file_sep)[[1]]
len_split = length(split_str) - 2
data_dir = paste(c(split_str[1:len_split],'data',''),collapse=file_sep)
# Check that the data directory contains the files for this weeks assignment
data_files = list.files(data_dir)
if(any(sort(data_files)!=sort(c('us_air_rev.csv','us_load_factor.csv')))){
cat("ERROR: DATA DIRECTORY NOT CORRECT\n")
cat(paste("data_dir variable set to: ",data_dir,collapse=""))
}
Problem 1 (Analytical Exercise)
Consider a simple AR(1) model:
Yt = αYt−1 + εt with εt ∼ N(0, σ
2
) for t = {1, . . . , T} and Y0 = 0
1. What is the distribution of Y1? What is the distribution of Y2?
2. What is the distribution of Yt for |α| < 1 as t → ∞.
3. What is the deffnition of stationarity? Explain why in this model we can check for stationarity by
looking at the mean and the variance of the Yt.
4. Suppose that α = 1. Why does this imply that the model is nonstationary? Can you think of a simple
transformation that makes the model stationary?
5. Now suppose that |α| < 1. Find a formula for the jth autocorrelation ρj = corr(Yt, Yt−j ).
26. Explain how we could use estimates of ρj for j = 1, 2, . . . to check whether some actual time series data
was generated by an AR(1) model like we one described above.
Problem 2 (Coding Exercise)
The problem will take you through a few tasks to familiarize yourself with R, as well as, some basic time
series concepts:
(a) Loading data into R
(b) Doing simple data analysis
(c) Doimg time series analysis
For this problem, we have pulled two seperate datasets from the FRED database, maintained by the Federal
Reserve Bank of Saint Louis (https://fred.stlouisfed.org/). The datasets cover the aggregate revenue
and load factor in domestic US ffights from 2000 to 2018. In the last two decades, airlines have begun
using sophisticated algorithms to increase capacity utilization of ffights (i.e. ffights tend to be more full).
Furthermore, during the same time period, airline revenues have increased. The point of this exercise will
be to understand the role of these productivity increases in “explaining” increased revenues in the airline
industry.
The two seperate datasets you will be working with are:
1. US Domestic Air Travel Revenue Passenger Mile (fflename = us_air_rev.csv) : this dataset contains
monthly data detailing the number of miles traveled by paying passengers in domestic US air travel.
2. US Domestic Air Travel Load Factor (fflename = us_load_factor.csv) : this dataset contains
monthly data detailing the percentage of seats fflled up (capacity utilitization) in domestic US air travel.
A Small Detour: Brief introduction to print statements We ask you to print a number of your
results in this exercise. In R, there are two different wants to print results:
1. print
2. cat
There are some deep programmatic differences underlying what each of these does, for our purposes we only
care about how easy to read your outputs are.
Printing Strings Let’s say you have a list of numbers, [4,5,6] and I want you to print out the following
statement:
The ffrst element of the list is: 4
The second element of the list is: 5
The third element of the list is: 6
Below I show you three ways to do so, the ffrst way simply uses print without any additional arguments.
The second way uses print with an additional argument, quote=False which gets rid of the quotes around
the strings. The third approaching, using cat, shows how this combines the second approach and has an
additional formatting feature that is useful for printing output.
## Define a list called x with 3 elements
x = c(4,5,6)
## Retrieve 1st, 2nd, 3rd element of list
first_elem = x[1] #1st element
second_elem = x[2] #2nd element
third_elem = x[3] #3rd element
## Create strings stating 'The XXXX element of the list is:'
first_str = 'The first element of the list is:'
3second_str = 'The second element of the list is:'
third_str = 'The third element of the list is:'
## Concatenate the list elements and the string to create the whole phrase
first_phrase = paste(first_str,first_elem,sep=' ')
second_phrase = paste(second_str,second_elem,sep=' ')
third_phrase = paste(third_str,third_elem,sep=' ')
## ~~ (1) Print without extra arguments ~~ ##
print('~~ (1) Print without extra arguments ~~')
print(first_phrase)
print(second_phrase)
print(third_phrase)
## ~~ (2) Print with extra argument turning off quotes ~~ ##
print('~~ (2) Print with extra argument turning off quotes ~~',quote=F)
print(first_phrase, quote=F)
print(second_phrase, quote=F)
print(third_phrase, quote=F)
## ~~ (3) Print without quotes and without trailing # ~~ ##
cat("\n")
cat("~~ (3) Print without quotes and without trailing # ~~\n")
cat(paste(first_phrase,"\n",sep=''))
cat(paste(second_phrase,"\n",sep=''))
cat(paste(third_phrase,"\n",sep=''))
[1] "~~ (1) Print without extra arguments ~~"
[1] "The first element of the list is: 4"
[1] "The second element of the list is: 5"
[1] "The third element of the list is: 6"
[1] ~~ (2) Print with extra argument turning off quotes ~~
[1] The first element of the list is: 4
[1] The second element of the list is: 5
[1] The third element of the list is: 6
~~ (3) Print without quotes and without trailing # ~~
The first element of the list is: 4
The second element of the list is: 5
The third element of the list is: 6
Printing Dataframes The main object you will be working with in R is called a dataframe (think an
excel spreadsheet). We will discuss more fully these objects in the following section. However, oftentimes you
will be asked to print out dataframes. In this case, using print is your best option. This is due to differences
between cat and print that are beyond the scope of this introduction.
(a) Loading Data
The ffrst thing we want you to do is to load both datasets: us_air_rev.csv and us_load_factor.csv into
R.
Please load data in the section below
## ~~ Problem 2: Part (a) Load Data into R ~~ ##
4## Change working directory to data directory
setwd(data_dir) # <- This is set in CODE SETUP section
# If you are having issues, you can manually
# set this variable to the folder where the data is
## Load Air Revenue Data ##
## Load Air Load Factor Data ##
There are two ways to view data that you have loaded into memory in R.
1. View only ffrst (or last few rows) using head (tails) commands
2. View the entire dataset in a seperate window using View commands
Note, for very large datasets it is not a good idea to use the View command as it is very memory (RAM)
intensive.
Other checks you always want to do when loading data includes:
1. Check the column names using colnames
2. Check the data types for each column using a loop and xxx
3. Check the dimension (number of rows and columns) using the dim command
We now want you to run the following checks on both of your loaded datasets:
(1) Print the column names.
(2) Print off the ffrst 20 rows.
(3) Print off the number of rows and columns.
(4) Print the data types of all the columns.
Note, for part (4) I have already built the for loop statement to get all the data types for each of the columns.
For those familiar with for loops in other environments, R has a built in set of apply functions that are
optimized for speciffc objects (lapply is optimized for lists, vapply is optimized for vectors etc). If you are
unfamiliar with for loops, give it a google.
## ~~ Problem 2: Part (a) Run Data Checks ~~ ##
(b) Doing simple data analysis
In the next part, we will have you doing some actual time series analysis. But generally we are interested in
decomposing time series into trend, seasonal and stochastic components. One clear form of seasonality is
month to month variation in the data. An “approximation” for trend components is to look at year to year
changes. We will have you investigate these below.
We now want you to do the following:
(1) Calculate the average revenue and load factor, by year. Do this two ways: (1) Using aggregate and
mean, (2) Using aggregate and sum.
(2) Calculate the average revenue and load factor, by month. Do this two ways: (1) Using aggregate and
mean, (2) Using aggregate and sum.
(3) Plot graphs for part (1) and (2) on the same plot, using your favorite plotting function. Note, you can
either use the built-in plot function or the popular external library ggplot2.
For parts (1) and (2), I want you to build a better understanding of using R. I am asking you to compute
averages using two different methods. In Method (1), you can use the built-in mean function to have R do
the work for you. In Method (2), you will do the average calculation yourself by summing over observations
and dividing by the number of observations.
## ~~ Problem 2: Part (b) Simple Data Analysis ~~ ##
5(c) Doing time series analysis
In R, there are already built-in functions that allow us to do these seasonality and trend decompositions
with much fewer lines of code. To do so, we must convert our data into time series objects. What seperates
a normal vector of data from a vector of time series data is that the latter has some time frequency of
observations. In our case, the time frequency is monthly.
We want now return to the main question of this section, how much does capacity utilizations explain increases
in airline revenue?
Fixing notation, we have:
• t ∈ {01/2000, . . . , 12/2018} = T is month-year combinations
• Revt is revenue for each month-year combination
• Loadt is load factor for each month-year combination
We now want you to do the following:
(1) Create a time series object using the ts command for each series. Be sure to specify the correct
frequency for the data.
(2) Plot an autocorrelation function between our two time series: {Revt}t∈T and {Loadt}t∈T .
(3) Run the following linear regression, reporting the coefficients and R2
:
Revt = α + βLoadt + ϵt
(4) Decompose both series into cyclical and trend components using the decompose command. Plot
seperately these cyclical and trend components for each of the series.
(5) Using the dataframe created in part (3), redo parts (2) and (3). What differs from part (2)? Why?
What can we conclude about the impact of capacity utilization changes on revenues?
## ~~ Problem 2: Part (b) Time Series Analysis ~~ ##
Problem 3 (Coding Exercise)
In class, you have learned about the Wold decomposition, a fundamental result in time series analysis. This
exercise will attempt to walk through Wold’s theorem in practice. We have provided simulated time series
data, in an .rda file called “ts_simulation.rda”, where Yt is the tth observation from our data. To open this
file, use the load command. The name of the dataframe is “sim”.
We now want you to do the following:
(1) Verify the stationarity of the process. Do this in two ways:
(a) “Heuristic” : show that the first-moment and second-moment do not depend on t.
(b) “Testing” : use a Dickey-Fuller test to test for stationarity. Interpret your results.
(2) Estimate three seperate autoregressive models: AR(1), AR(3) and AR(6). For each of the seperate
models, retrieve the residuals, ϵˆ{t,p}, where p is the order of the AR process. Using each set of residuals
of the AR process, estimate an MA(2) model, where ηˆ{t,p,q} are the residuals of this second step.
Verify whether the assumptions of Wold are violated:
Corr[ˆϵ{t,p}Ys] = 0 such that s < t
E[η{t,p,q}] = 0
V ar(η{t,p,q}) = σ
2
(3) To find the right ARMA(p,q) process, we add new lags (increase p), estimate our model, use an
information criteria to determine the increase in fit and stop once new models do not improve fit. To
simplify the problem, assume q = 2. Build a series of ARMA(p,q) models, using the Akaike Information
Criteria (AIC) to find the right p. (Note: A for loop over p would be a good idea).
6## ~~ Problem 3: Wold-Decomposition ~~ ##
7