Wrapping up Graphics and Introducing Time Series
STAT 4360
Friday, April 19
I) Multiple Plots + Window Settings
You can change the number of plots in the graphing window using par(mfrow=c(#rows,#columns)). The par() function can also rotate the axis labels and increase axis margins. The layout() function allows for more window options and is more flexible than par().
layout(matrix(c(1,3,2,4),2,2,byrow=T))In this piece of code, we specify four quadrants (two rows and two columns) of equal size for the plotting window. The numeric vector c(1,3,2,4) specifies the order in which the plot is created. If you put a 0 inside of this vector, this quadrant will be left empty. Aditional function inputs can change the widths of the rows (heights=c(1,2)) and the widths of the columns (widths=c(3,1)).
## Changing plotting window
layout(matrix(c(1,3,2,4),2,2,byrow=T))
## Filling with plots
plot(1:10,1:10)
plot(1:10+10,11:20)
plot(1:10+20,21:30)
plot(1:10+30,31:40)
In the following example, we change with width of the columns.
## Changing plotting window
layout(matrix(c(1,3,2,4),2,2,byrow=T),widths=c(3,1))
## Filling with plots
plot(1:10,1:10)
plot(1:10+10,11:20)
plot(1:10+20,21:30)
plot(1:10+30,31:40)
In the following example, we change the width of the rows.
## Changing plotting window
layout(matrix(c(1,3,2,4),2,2,byrow=T),heights=c(1,2))
## Filling with plots
plot(1:10,1:10)
plot(1:10+10,11:20)
plot(1:10+20,21:30)
plot(1:10+30,31:40)
You can also change both at once.
## Changing plotting window
layout(matrix(c(1,3,2,4),2,2,byrow=T),widths=c(3,1),heights=c(1,2))
## Filling with plots
plot(1:10,1:10)
plot(1:10+10,11:20)
plot(1:10+20,21:30)
plot(1:10+30,31:40)
If you would like to clear and return to default plot settings, run def.off(), or device off.
J) Saving Graphics
- From the graphical user interface (GUI): Export -> Save as Image… and then select the dexired format.
- You can also save graphics to the working directory (or another specified directory) using functions. For example,
pdf('mygraph.pdf')tells R you’re dumping the following graphics into a PDF file. Similar functions include:win.metafile(mygraph.wmf),png(mygraph.png),jpeg(mygraph.jpeg),bmp(mygraph.bmp), andpostscript(mygraph.ps)perform similar roles.
ASSIGNMENT 12!!
Statistics in R
Based on the pre-reqs for this course, everyone should already know how to do simple regression and ANOVA in R. These will be small crash courses in statistics that you might not have seen.
Time Series
- A time series is a sequence of observations taken over (equally-spaced) time
- Goal: predict future observations given previous observations
- Models given here assume constant mean and independent N(0,\(\sigma^2\)) errors.
Some steps for working with time series:
- Plot the series. Are \(\mu\) and \(\sigma^2\) constant?
- Examine autocorrelation function (ACF) and partial autocorrelation function (PACF) plots. We’ll talk about these more later.
- Fit a time series model.
- Forecast with the model.
Four potential time series models:
- White Noise = uncorrelated. Recall that independent data are uncorrelated, but uncorrelated data are not necessarily independent.
Model: \(x_t=\mu+\epsilon_t\)
## White Noise
obs=100
set.seed(3)
x=rnorm(obs)
## This code will create a time series object
ts.x=ts(x)
## Now, we plot it + ACF + PACF
layout(matrix(c(1,1,2,3),2,2,byrow=T))
plot(ts.x,type='b',main='White Noise Series',pch=20,ylab=expression(x(t)))
acf(ts.x,main='White Noise ACF',ylim=c(-1,1))
pacf(ts.x,main='White Noise PACF',ylim=c(-1,1))
Things to note:
- Nothing special in the series itself. Look at it through the window and you will see the same/similar behavoir everywhere (“stationarity”).
- ACF: correlation of the series with itself at different lags. The ACF = 1 at lag 0 indicating that you’re perfectly correlated with yourself. \(cor(x_t,x_{t+h})\), where \(h\) indicates the lag. The dashed lines represent 95% confidence interval at each lag. Note that there is no significant correlation at \(h>0\).
- PACF: coefficient when regressing on one more \(x_t\). This will make a lot more sense in a minute! The PACF plot does not show lag 0. Dashed lines are 95% confidence interval at lags. Again, there is no significance at any \(h\).
- Autoregressive = regress on yourself. Assume \(\mu=0\). Or, substitute \(x_t-\mu\) for all \(t\) to essentially center the time series.
In an \(AR(1)\) model, we include 1 lag (i.e. today depends on yesterday). The model form: \(x_t=\phi x_{t-1} + \epsilon_t\).
In an \(AR(p)\) model, we include \(p\) lags (i.e. today depends on the previous \(p\) days). The model form \(x_t=\phi_1 x_{t-1}+\dots+\phi_p x_{t-p}+\epsilon_t\). The \(\phi\) coefficients are the PACF at different lags. “Regress on one more - what’s the coefficient?” The \(\phi\) coefficients must be between 0 and 1. If they are all zero, then you have a white noise series and vice versa.
## AR(1)
burn=400
phi=0.5 ## Play around with this parameter
set.seed(3)
## Generate sequence
x=numeric(burn+obs) ## Empty vector
x[1] = rnorm(1)
for(i in 2:(obs+burn)){
x[i] = phi*x[i-1]+rnorm(1)
}
x=x[(burn+1):(burn+obs)]
## Creating a time series object
ts.x=ts(x)
## Plot
layout(matrix(c(1,1,2,3),2,2,byrow=T))
plot(ts.x,type='b',main='AR(1) Series',pch=20,ylab=expression(x(t)))
acf(ts.x,main='White Noise ACF',ylim=c(-1,1))
pacf(ts.x,main='White Noise PACF',ylim=c(-1,1))
Some things to note:
- If ACF decays exponentially with lag this implies you have an \(AR(p)\) model.
- PACF at lag 1 is about equal to 0.5, all others about equal to 0. These values should decay, too: Wednesday should be more strongly related to today than Tuesday.
## Fitting the AR(1) model
ar1 = arima(ts.x,order=c(1,0,0))
print(ar1)##
## Call:
## arima(x = ts.x, order = c(1, 0, 0))
##
## Coefficients:
## ar1 intercept
## 0.4735 0.3027
## s.e. 0.0888 0.2096
##
## sigma^2 estimated as 1.24: log likelihood = -152.76, aic = 311.52