# The first time you use R (and perhaps after that as well) you will want to enter these three commands:
install.packages("astsa")
require(astsa)
help(package = astsa)




# This is a brief introduction to R
# Open R, and then work through these commands

#   '#'  = comment
#   '<-' or '=' = affect; defines the left hand side to be the right hand side
#   '>'  = prompt of R

# basic operations - enter these; see what you get: 

7+5

7-5

7*5

7/5

7/0

0/0

# use the function c() for creating a vector 
x = c(1,2,3) 
x  # response: [1] 1 2 3
x = 1:3 # Gives the same thing
x

x = c(1,2,3) + c(2,3,4) 
x # response: [1] 3 5 7

c(1,2,3)*c(2,3,4) # response: [1]  2  6 12

c(1,2,3)*2 # response: [1] 2 4 6

c(1,2,3)/2 # response: [1] 0.5 1.0 1.5

# create a vector, called 'height', containing 10 values
height <- c(72,70,71,66,73,66,72,72,74,75)
height # response: [1] 72 70 71 66 73 66 72 72 74 75

mean(height)              # compute sample average

median(height)            # compute sample median

mean(height, trim=0.1)    # trimmed mean.

var(height)               # variance
sum((height - mean(height))^2)/(length(height)-1)
				# should give the same thing - why?

sqrt(var(height))         # standard deviation

height[2]                 # the second element of the vector

       
# generate a vector 
rep(2,3)               # a vector of 2's of length 3: [1] 2 2 2

seq(1,10, length=21)   # a vector of length 21; equally spaced entries spanning [1,10]

# Two time series already in R are 'sales' and 'lead' - one measures sales and the other 
# is allegedly a 'leading indicator' that gives advance information about sales.
x = lead
  # Is it already a time series?
is.ts(x)
y = sales
sales  #Print out the data in 'sales'
SalesAndLead = cbind(sales, lead) # Create a matrix containing both
SalesAndLead

plot(sales)
plot(lead)
plot(x=lead, y=sales)
# This plot will be UGLY, because "lead" and "sales" are time series, and
# a plot of one against the other will have the points connected by lines.
# We really wanted to plot them as one vector of numbers against the other:
plot(x=as.vector(lead), y=as.vector(sales))
# Even nicer:
plot(x=as.vector(lead), y=as.vector(sales), xlab = "lead", ylab = "sales")


acf(sales)  # Estimate of the atocorrelation function
acf(lead, 36)  #What did '36' do?

win.graph() # Opens a new graphics window without closing the current one
par(mfrow=c(2,1)) # 2 by 1 panel of plots
ts.plot(sales, lead, lty=2:1, col=1:2, lwd=1:2, main="sales (top) and leading indicator")
ts.plot(sales, lead, lty=1:1, col=2:3, lwd=1:2, main="sales (top) and leading indicator")  # What has changed?
help(par)  # help on commands like 'lty' (= line type) and 'lwd' (= line width)

win.graph()
par(mfrow=c(2,1))
fit = lm(sales ~ lead)  # A simple linear regression of sales on the leading indicator
summary(fit)
resids = fit$resid 
ccf(lead,resids)  # Seems to say that the earlier values of 'lead' will help explain the remaining changes in sales
			# It wouldn't be called a 'leading' economic indicator otherwise, would it?
ts.plot(resids)  # Much remains to be explained

## For this next step, first click on "Packages" in the menu at the top of the R screen.
## Choose the Canadian mirror, and then when prompted ask to have the 'dynlm' package installed.
## If you've done this once already you shouldn't have to do it again (on the same computer).

library(dynlm)

# This function dynlm fits linear models with lagged series; 
# it addresses the problem that the lagged series will have differing lengths

## Regress Y = Sales on X = Lead, at various lags.
par(mfrow=c(2,2)) # 2 by 2 panel of plots, by rows
ts.plot(resids) # Repeat the previous two plots
abline(h=0) # What does this do?
ccf(lead,resids)
  # Regress Sales on past (3 - 9 months, chosen by trial and error) leads
fit2 = dynlm(sales ~ L(lead,3:9)) # The same as sales ~ L(lead,3) + ... + L(lead,9))
resids2 = fit2$resid
ts.plot(resids2)
abline(h=0)
ccf(lead,resids2) # Residuals are much smaller; no significant cross-correlations
summary(fit2)



# Using 'scripts':  What you are reading right now is an example of a 'script' - a list of commands
# to be carried out by R.  In your work you will typically have a script for each assigned question, etc.
# To open a blank script, click on File -> New script in the R menu at the top of the screen.
# To save the script for future use, click on the 'save' icon and follow the instructions.  
# Next time you want to use it, double click on it - R should open in that same directory (if not, 
# you might have to hunt around for the script, after clicking on File -> Open script.)
# To run the command in the script, save it and then click on Edit -> Run all.  
# (Try it - to run THIS script, save it to a directory somewhere, double-click on it to re-open R from 
# that directory, then open it and run it.)
# Alternatively, one will often just copy and paste a few commands, from the script, 
# into the Command window ("R console") in order to try them out.