### CLASS EEG #######

## BEFORE PROCEEDING DO THE FF
## Load ts package
## Change directory to where the datasets are located

## Columns correspond to EEGS at the following locations C3 P3 F7 T3 T5
eegA = read.table("EEGApre.txt");
eegB = read.table("EEGBpost.txt");
dim(eegA)
dim(eegB)

T3pre=eegA[, 4];
T3post=eegB[,4];

## compare time plots of pre and during/post
par(mfrow=c(3,1))
plot(T3pre, type="l"); title(paste("Pre seizure"));
plot(T3pre, type="l", ylim=c(-250, 450)); title(paste("Pre seizure"));
plot(T3post, type="l", ylim=c(-250, 450)); title(paste("During/Post seizure"));


par(mfrow=c(3,1));
plot(T3pre, type="l");
acf(as.ts(T3pre), lag.max=100);
acf(as.ts(T3pre), type="partial", lag.max=100);

## AR(2) seems to be appropriate. Let's try AR(2), AR(3), AR(4) and AR(5) and compare

x=as.ts(T3pre - mean(T3pre));
Mod2 = arima(x, order=c(2,0,0), include.mean=F);
Mod3 = arima(x, order=c(3,0,0), include.mean=F);
Mod4 = arima(x, order=c(4,0,0), include.mean=F);
Mod5 = arima(x, order=c(5,0,0), include.mean=F);
Mod6 = arima(x, order=c(6,0,0), include.mean=F);

## Function: for a given model, compute aic, aic, bic
infocrit = function(modout){

aic = log(modout$sigma2) + (2*(length(modout$coef)+1))/length(modout$resid);
aicc = log(modout$sigma2) + (2*(length(modout$coef)+1))/(length(modout$resid) - 
(length(modout$resid) + 2))
bic = log(modout$sigma2) + (length(modout$coef)+1)*log(length(modout$resid))/length(modout$resid);

return(c(aic, aicc, bic))}


## Model2, 3, 4,5,6:
infocrit(Mod2)
infocrit(Mod3)
infocrit(Mod4)
infocrit(Mod5)
infocrit(Mod6)
 
## According to BIC: AR(5)

## Do the likelihood ratio test if AR(6) can be simplified into AR(5) 
## Ho: phi(6) = 0 vs H1: phi(6) ne 0
## Test Statistic
chisqstat = 2*(Mod6$loglik - Mod5$loglik)
 ## Cut-off value
qchisq(0.95, 1)
## Conclusion: Reduce AR(6) to AR(5)

## Can we reduce AR(5) to AR(4)?
## Ho: phi(5) = 0 vs H1: phi(5) ne 0
chisqstat = 2*(Mod5$loglik - Mod4$loglik)
 ## Cut-off value
qchisq(0.95, 1)
## Conclusion: Don't reduce. Select AR(5) as the final model

## Check residuals for AR(5)
par(mfrow=c(3,1))
plot(Mod5$residuals, type="l");
acf(as.ts(Mod5$residuals), lag.max=200);
acf(as.ts(Mod5$residuals), type="partial", lag.max=200);
## Conclusion: AR(5) is a reasonable model. 

## Check the Box-Lyung tests for independence 
K = 27; 
acf.Mod5 = acf(as.ts(Mod5$residuals));
QLB = 512*(514)*sum(((acf.Mod5$acf[2:28])^2)/(512-acf.Mod5$lag[2:28]));
pval = 1 - pchisq(QLB, df=27)
## No evidence against independence

## Would AR(2) have been a good model?
par(mfrow=c(3,1))
plot(Mod2$residuals, type="l");
acf(as.ts(Mod2$residuals), lag.max=200);
acf(as.ts(Mod2$residuals), type="partial", lag.max=200);
## Conclusion: Not quite!. 

## AR(5) reasonable model 
## Examine the coefficients and their standard errors
Mod5$coef
sqrt(diag(Mod5$var.coef))

Mod5$coef/sqrt(diag(Mod5$var.coef))

## Go back to the time plot
par(mfrow=c(2,1))
plot(T3pre, type="l", ylim=c(-200,400))
title(paste("EEG T3 channel Pre Seizure"));
plot(T3post, type="l")
title(paste("EEG T3 channel Post Seizure"));

########################################################################
########################################################################
########################################################################
## Investigate post seizure correlation properties
##
par(mfrow=c(3,1));
plot(T3post, type="l");
acf(as.ts(T3post), lag.max=100);
acf(as.ts(T3post), type="partial", lag.max=100);


x=as.ts(T3post - mean(T3post));
Mod22 = arima(x, order=c(2,0,0), include.mean=F);
Mod32 = arima(x, order=c(3,0,0), include.mean=F);
Mod42 = arima(x, order=c(4,0,0), include.mean=F);
Mod52 = arima(x, order=c(5,0,0), include.mean=F);
Mod62 = arima(x, order=c(6,0,0), include.mean=F);

infocrit(Mod22)
infocrit(Mod32)
infocrit(Mod42)
infocrit(Mod52)
infocrit(Mod62)

## BIC also chooses AR(5)

## Check residuals for AR(5)
par(mfrow=c(3,1))
plot(Mod52$residuals, type="l");
acf(as.ts(Mod52$residuals), lag.max=200);
acf(as.ts(Mod52$residuals), type="partial", lag.max=200);
## Conclusion: AR(5) is a reasonable model. 

## Check the Box-Lyung tests for independence 
K = 27; 
acf.Mod5 = acf(as.ts(Mod52$residuals));
QLB = 512*(514)*sum(((acf.Mod5$acf[2:28])^2)/(512-acf.Mod5$lag[2:28]));
pval = 1 - pchisq(QLB, df=27)
## No evidence against independence
##  Quite contradictory to the acf

## Show the estimates of 2 models

Mod5
Mod52

## Check the roots of the AR polynomial functions
polyroot(c(1, -Mod5$coef))
polyroot(c(1, -Mod52$coef))

## Check the magnitudes of the roots
abs(polyroot(c(1, -Mod5$coef)))
abs(polyroot(c(1, -Mod52$coef)))


################## Spectral Analysis #############################
################## Spectral content of the signals ##################

fft1 = fft(T3pre);
period1 = ((abs(fft1))^2)/length(T3pre);

fft2 = fft(T3post);
period2 = ((abs(fft2))^2)/length(T3pre);

par(mfrow=c(2,1));
plot(period1[1:256], type="l", lwd=3);
title(paste("Pre-seizure"));
plot(period2[1:256], type="l", lwd=3);
title(paste("During/Post seizure"));

par(mfrow=c(1,1));
ymax=max(max(period1, period2));
plot(period1[1:256], type="l", col=1,ylim=c(0, ymax), lwd=3);
lines(period2[1:256], col=2, lwd=3)
title(paste("Periodograms Pre vs Post"));


par(mfrow=c(2,1));
ymax=max(max(period1, period2));
plot(period1[1:256], type="l", col=1,ylim=c(0, ymax), lwd=3);
lines(period2[1:256], col=2, lwd=3)
title(paste("Periodograms Pre vs Post"));
normperiod1 = 100*period1/sum(period1);
normperiod2 = 100*period2/sum(period2);
plot(normperiod1[1:256], type="l", lwd=3);
lines(normperiod2[1:256], col="red", lwd=3)
title(paste("Normalized Periodograms Pre vs Post"));