Stats 429 - Time Series Analysis

Intro to R

--------------------------------------

Outline

I. Basics

II. General Commands (graphs, reading data sets, etc.)

III. Create simple functions

IV. Time series functions

----------------------------------------

I. BASICS

1. Always create a file using any editor to help save time by not 
having to re-type commands.

2. Help - use pull down menu or type in the command line 

> help.search("whateveritisyouneed")

3. Use # to comment out a line.

4. Remember to "load" the ts package.

------------------------------------------

II. GENERAL COMMANDS

- load a data set
- vector and matrix manipulation
- plots

## Change directory to where the datasets are located
## Columns correspond to EEGS at the following locations C3 P3 F7 T3 T5

eegpre = read.table("EEGApre.txt");
eegpost = read.table("EEGBpost.txt");

#dim(eegpre)
#dim(eegpost)

C3pre=eegpre[, 1];
C3post=eegpost[,1];

P3pre=eegpre[, 2];
P3post=eegpost[,2];

T3pre=eegpre[, 4];
T3post=eegpost[,4];

## Compare time plots of pre and during/post
## Plots (3 rows and 2 columns)

par(mfrow=c(3,2));

## Plot 1
Tlen = length(C3pre);
xt = seq(1/Tlen,1, by=(1/Tlen));
y1t = C3pre; 
y2t = C3post;
plot(xt, y1t, xlab="Time", ylab="EEG", ylim=c(-300, 500)); title(paste("C3 Pre seizure"));
plot(xt, y2t, xlab="Time", ylab="EEG", ylim=c(-300, 500)); title(paste("C3 Post seizure"));

## Plot 2
plot(P3pre, type="l", lwd=2); title(paste("P3 Pre seizure"));
plot(P3post, type="l", lwd=3); title(paste("P3 During/Post seizure"));

## Plot 3
plot(T3pre, type="l"); title(paste("T3 Pre seizure"));
plot(T3post, type="l"); title(paste("T3 During/Post seizure"));

## To save the plots in one file

- click on the graphics window (activate it)
- click on right button of mouse, save as postscript

--------------------------------------------

III. Create simple functions

## 1. Create a function that computes the correlation coefficient between 
## vectors Y and x

correst<-function(y,x)
{
## Input: y and x
## Output: corrout

xnew = x - mean(x);
ynew = y - mean(y);
sumxy = sum(xnew*ynew);
varxy = sumxy/(length(y) -1);

varx = var(x);
vary = var(y);

corrout = varxy/sqrt(varx*vary);

## Output is corrout 
corrout}

## Example
x = rnorm(100, 0, 1);
y = 2*x + rnorm(100, 0, 1);
par(mfrow=c(1,1)); plot(x, y, xlab="x", ylab="y"); title(paste("Y vs. X"));
mycorrest = correst(x,y);
mycorrest


## 2. Create a function that can generate AR(p) time series

arts = function(T, arparm, noisesd)
{
## Input: T = length of time series
## 	  arparm = coefficients of the AR
##	noisesd = sdt dev of the noise
## Output: tsout

p = length(arparm);
templength = T+1000+p;
innov = rnorm(templength, mean=0, sd=noisesd);

tsout = rep(0, templength);
tsout[1:p] = innov[1:p];
for (k in (p+1):templength){
temp = tsout[(k-p):(k-1)];
tsout[k] = sum(arparm*temp) + innov[k];
}
tsout = tsout[(templength-T+1):templength]
tsout} 

## Example
tsdataset1 = arts(T=500, arparm = c(0.6, -0.3), noisesd=1);
tsdataset2 = arts(T=500, arparm = c(0.9), noisesd=1);
tsdataset3 = arts(T=500, arparm = c(-0.9), noisesd=1);

## Plot the time series

par(mfrow=c(3,1));
plot(tsdataset1, type="l"); title(paste("AR(2) (0.6, -0.3)"));
plot(tsdataset2, type="l"); title(paste("AR(1) (0.9)"));
plot(tsdataset3, type="l"); title(paste("AR(1) (-0.9)"));
----------------------------------------

IV. Some Time Series Functions

1. arima
2. acf
3. diff
4. filter
5. fft