## Correlations and Models 

##### White Noise Process on Uniform(-a,a)
T=100; a=3;
x=runif(T, (-a), a); 

par(mfrow=c(2,1));
hist(x);
ts.plot(x);

##### Gaussian White Noise Process
T=100; mu=0; sig2 = 4; 
x=rnorm(T, mu, (sqrt(sig2)));

par(mfrow=c(2,1));
hist(x);
ts.plot(x);

## Correlation between ynew=X_{t} and xnew=X_{t-1}
ynew = x[2:T];
xnew = x[1:(T-1)];
par(mfrow=c(1,1));
plot(xnew, ynew, type="p");
cor(xnew, ynew)

##### First Order Moving Average MA(1) with Guassian white noise
T=100; mu=0; sig2=1; theta=1;
x=c(1:101);
wn=rnorm((101), mu, (sqrt(sig2)));
for (k in 2:101){
x[k] = wn[k] + theta*wn[(k-1)];
}
x=x[2:101];

##### Built-in R function for generating MA time series data
x = arima.sim(T, model=list(ma=c(theta)));


par(mfrow=c(1,1)); 
ts.plot(x); 

par(mfrow=c(2,2)); 
xnew=x[1:(T-1)]; ynew=x[2:T]; plot(xnew, ynew, type="p"); title(paste("Lag 1"));
cov(xnew, ynew)
cor(xnew, ynew)

xnew=x[1:(T-2)]; ynew=x[3:T]; plot(xnew, ynew, type="p"); title(paste("Lag 2"));
cov(xnew, ynew)
cor(xnew, ynew)  

xnew=x[1:(T-3)]; ynew=x[4:T]; plot(xnew, ynew, type="p"); title(paste("Lag 3"));
cov(xnew, ynew)
cor(xnew, ynew) 

xnew=x[1:(T-4)]; ynew=x[5:T]; plot(xnew, ynew, type="p"); title(paste("Lag 4"));
cov(xnew, ynew)
cor(xnew, ynew) 

## The sample auto-correlation 
par(mfrow=c(1,1));
acf(x);

##### First Order Auto-Regressive Model AR(1) with Gaussian white noise

T=100; burnin=100; mu=0; sig2=1; phi=0.9;
x=c(1:(T+burnin));
wn=rnorm((T+burnin), mu, (sqrt(sig2))); 
x[1]=0;
for (k in 2:(T+burnin)){
x[k] = phi*x[(k-1)] + wn[k];}

plot(x, type="l");
x = x[(burnin+1):(T+burnin)];

plot(x, type="l"); 

##### R built in function for generating AR time series data 
x = arima.sim(T, model=list(ar=c(phi)));

par(mfrow=c(2,2)); 
xnew=x[1:(T-1)]; ynew=x[2:T]; plot(xnew, ynew, type="p"); title(paste("Lag 1"));
cov(xnew, ynew)
cor(xnew, ynew)

xnew=x[1:(T-2)]; ynew=x[3:T]; plot(xnew, ynew, type="p"); title(paste("Lag 2"));
cov(xnew, ynew)
cor(xnew, ynew)  

xnew=x[1:(T-3)]; ynew=x[4:T]; plot(xnew, ynew, type="p"); title(paste("Lag 3"));
cov(xnew, ynew)
cor(xnew, ynew) 

xnew=x[1:(T-4)]; ynew=x[5:T]; plot(xnew, ynew, type="p"); title(paste("Lag 4"));
cov(xnew, ynew)
cor(xnew, ynew) 

## The sample auto-correlation 
par(mfrow=c(1,1));
acf(x);

## Note the difference between plots for phi=0.9 vs phi=-0.9
T=1024; burnin=100; mu=0; sig2=1; phi=0.9;
x=c(1:(T+burnin));
wn=rnorm((T+burnin), mu, (sqrt(sig2))); 
x[1]=0;
for (k in 2:(T+burnin)){
x[k] = phi*x[(k-1)] + wn[k];}
x = x[(burnin+1):(T+burnin)];
par(mfrow=c(2,1));
plot(x, type="l"); title(paste("AR(1), phi=+0.9"));


T=1024; burnin=100; mu=0; sig2=1; phi=-0.9;
x=c(1:(T+burnin));
wn=rnorm((T+burnin), mu, (sqrt(sig2))); 
x[1]=0;
for (k in 2:(T+burnin)){
x[k] = phi*x[(k-1)] + wn[k];}
x = x[(burnin+1):(T+burnin)];
plot(x, type="l"); title(paste("AR(1), phi=-0.9"));


######### BUILT In R functions

### generate N=100 from ARMA(1,1) with parameters 0.5 and 0.9: 
N=100;
tsdata = arima.sim(N, model=list(ar=c(0.5), ma=c(0.9)));
plot(tsdata);
 
### generate N=100 from AR(1) 0.9
N=100;
tsdata1 = arima.sim(N, model=list(ar=c(0.5)));
plot(tsdata1);

### Estimate the ACF
acf(tsdata1, plot=T)

### Estimate the partial ACF
pacf(tsdata1,plot=T)
 

### generate N=100 from MA(1) 0.9
N=100;
tsdata2 = arima.sim(N, model=list(ma=c(0.9)));
plot(tsdata2);

### Estimate the ACF
acf(tsdata2, plot=T)

### Estimate the partial ACF
pacf(tsdata2,plot=T)

### Plot them side by side
par(mfrow=c(3,2))
plot(tsdata1);
plot(tsdata2);
acf(tsdata1, plot=T)
acf(tsdata2, plot=T)
pacf(tsdata1, plot=T) 
pacf(tsdata2, plot=T)