## Spectral Analysis Fall 2006


source("C:/Documents and Settings/ombao/My Documents/CLASS/Fall 2006/Stats429-Fall06/R Files/armaspec.R")
source("C:/Documents and Settings/ombao/My Documents/CLASS/Fall 2006/Stats429-Fall06/R Files/periodsmoother.R")
source("C:/Documents and Settings/ombao/My Documents/CLASS/Fall 2006/Stats429-Fall06/R Files/powertf.R")


########## Linear Regression #######################################

## Time series = a mixture of components of some trigonometric basis
## Ingredients = The sine waves, cosine waves, Fourier complex exponentials           
## Random coefficients = amount or proportion of each of the waveform


## Fourier basis
T = 256;
timet = c(0:(T-1));
freqindex = c(-(T/2-1):(T/2));
FB = matrix(ncol=T,nrow=T);
for(k in 1:T){
tempfreq = 2*pi*freqindex[k]/T;
FB[,k] = complex(real=cos(tempfreq*timet), imag=sin(tempfreq*timet));
}

## The "Low Frequency Waveforms" 
tempfreq = 2*pi*2/T;
fwave = complex(real=cos(tempfreq*timet), imag=sin(tempfreq*timet));
par(mfrow=c(1,1));
plot((1:T), Re(fwave), type="l", ylim=c(-1.1,1.1), col=1);
lines((1:T), Im(fwave), type="l", col=2);

tempfreq = 2*pi*-1/T;
fwave = complex(real=cos(tempfreq*timet), imag=sin(tempfreq*timet));
par(mfrow=c(1,1));
plot((1:T), Re(fwave), type="l", ylim=c(-1.1,1.1), col=1);
lines((1:T), Im(fwave), type="l", col=2);


## The "Middle Frequency Waveforms"
tempfreq = 2*pi*20/T;
fwave = complex(real=cos(tempfreq*timet), imag=sin(tempfreq*timet));
par(mfrow=c(1,1));
plot((1:T), Re(fwave), type="l", ylim=c(-1.1,1.1), col=1);
lines((1:T), Im(fwave), type="l", col=2);

## The "High Frequency Waveforms"
tempfreq = 2*pi*20/T;
fwave = complex(real=cos(tempfreq*timet), imag=sin(tempfreq*timet));
par(mfrow=c(1,1));
plot((1:T), Re(fwave), type="l", ylim=c(-1.1,1.1), col=1);
lines((1:T), Im(fwave), type="l", col=2);


## Are the Fourier waveforms orthogonal? 
## Check if (FB)' (FB) is a diagonal matrix
## Matrix whose (j,k) element is the inner product between 
## j-th wave and k-th wave.
IPmatrix = t(Conj(FB))%*%FB;

## Then examine the off diagonals of IP. One way is to plot 
## each row. Check if only the r-th element of row r is 
## large and the rest are essentially zero.

## Normalize the Fourier basis 
FB = (1/sqrt(T))*FB;


########################################
## Mixing the ingredients

## At frequency k=0, the cosine wave has a coefficient that 
## is randomly generated from a normal distribution with mean 0
## and variance v(0).

## At frequency k=1, ..., T/2 -1, each of the cosine and sine 
## waves has coefficient from normal with mean 0 and variance 
## v(k)/2. We say that at frequency k=1, ..., N/2-1, the total 
## variance for sine and cosine part is v(k). 

## At frequency k=T/2, the cosine wave has a coefficient that 
## is randomly generated from a normal distribution with mean 0
## and variance v(N/2).

## All the coefficients are independent (though not necessarily 
## identically distributed since the variance may depend on 
## the frequency index k).

############ Example 1 ###################
## A mixing with predominantly low frequencies
## Power is concentrated on the lowest 20 of frequencies
## The variance function is v(k)

## The mixing variance v(k)
V = c(1:(T/2 + 1));
V[1:20] = 100; 
V[21:50] = 30;
V[51:(T/2+1)] = 1;
par(mfrow=c(1,1));
plot(V, type="l", xlab="Frequency index", ylab="Spectrum");

par(mfrow=c(1,1))
vmat = matrix(V);
image(t(-vmat), axes=FALSE, col=heat.colors(12), xlab="Time", ylab="Frequency");
title(paste("Time-Frequency Plot of Variance of process"));


########## Generate the Random Coefficients
RandCoeff = c(1:T); 

PosCoeff = c(1:(T/2+1));
PosCoeff[1] = rnorm(1,0, sqrt(V[1]));
PosCoeff[(T/2+1)] = rnorm(1,0, sqrt(V[T/2+1]));
for(k in 2:(T/2)){
Realpart = rnorm(1,0, sqrt(V[k]/2)); 
Imagpart = rnorm(1,0, sqrt(V[k]/2));
PosCoeff[k] = complex(real=Realpart, imag=Imagpart);
}

RandCoeff[(T/2):T] = PosCoeff;
RandCoeff[1:(T/2-1)] = Conj(PosCoeff[(T/2):-1:2]);

### plot the real part of the coefficents
par(mfrow=c(1,2)); 
plot(Re(RandCoeff), type="l");
plot(Im(RandCoeff), type="l");

### plot the squared coefficients
par(mfrow=c(1,1)); 
Pgram = (abs(RandCoeff))^2;
plot(Pgram, type="l")

### Compare this with the variance function 
par(mfrow=c(1,1)); 
ymax = max(Pgram, V); ymin = min(Pgram, V); 
plot((1:(T/2+1)), V[1:(T/2+1)], ylim=c(ymin, ymax), xlab="Frequency", ylab="Var", col=1); 
lines(Pgram[(T/2):T], col=2)

### The time series 
tsdata = FB%*%RandCoeff;
plot(Im(tsdata))
tsdata = Re(tsdata);
par(mfrow=c(1,1)); 
plot(tsdata, type="l")


############ Example 2 ###################
## A mixing with predominantly high frequencies
## Power is concentrated on the highest 20 of frequencies
## The variance function is v(k)

## The mixing variance v(k)
V = c(1:(T/2 + 1));
V[1:29] = 1; 
V[30:109] = 30;
V[110:129] = 100;
par(mfrow=c(1,1));
plot(V, type="l", xlab="Frequency index", ylab="Spectrum");

par(mfrow=c(1,1))
vmat = matrix(V);
image(t(-vmat), axes=FALSE, col=heat.colors(12), xlab="Time", ylab="Frequency");
title(paste("Time-Frequency Plot of Variance of process"));


########## Generate the Random Coefficients
RandCoeff = c(1:T); 

PosCoeff = c(1:(T/2+1));
PosCoeff[1] = rnorm(1,0, sqrt(V[1]));
PosCoeff[(T/2+1)] = rnorm(1,0, sqrt(V[T/2+1]));
for(k in 2:(T/2)){
Realpart = rnorm(1,0, sqrt(V[k]/2)); 
Imagpart = rnorm(1,0, sqrt(V[k]/2));
PosCoeff[k] = complex(real=Realpart, imag=Imagpart);
}

RandCoeff[(T/2):T] = PosCoeff;
RandCoeff[1:(T/2-1)] = Conj(PosCoeff[(T/2):-1:2]);


### plot the real part of the coefficents
par(mfrow=c(1,2)); 
plot(Re(RandCoeff), type="l");
plot(Im(RandCoeff), type="l");

### plot the squared coefficients
par(mfrow=c(1,1)); 
Pgram = (abs(RandCoeff))^2;
plot(Pgram, type="l")

### Compare this with the variance function 
par(mfrow=c(1,1)); 
ymax = max(Pgram, V); ymin = min(Pgram, V); 
plot((1:(T/2+1)), V[1:(T/2+1)], ylim=c(ymin, ymax), xlab="Frequency", ylab="Var", col=1); 
lines(Pgram[(T/2):T], col=2)

### The time series 
tsdata = FB%*%RandCoeff;
plot(Im(tsdata))
tsdata = Re(tsdata);
par(mfrow=c(1,1)); 
plot(tsdata, type="l")



############ Example 3 ###################
## All frequencies of the same weight 
##

V = c(1:(T/2 + 1));
V[1:(T/2+1)] = 1;
par(mfrow=c(1,1));
plot(V, type="l", xlab="Frequency index", ylab="Spectrum");

par(mfrow=c(1,1))
vmat = matrix(V);
image(t(-vmat), axes=FALSE, col=heat.colors(12), xlab="Time", ylab="Frequency");
title(paste("Time-Frequency Plot of Variance of process"));


########## Generate the Random Coefficients
RandCoeff = c(1:T); 

PosCoeff = c(1:(T/2+1));
PosCoeff[1] = rnorm(1,0, sqrt(V[1]));
PosCoeff[(T/2+1)] = rnorm(1,0, sqrt(V[T/2+1]));
for(k in 2:(T/2)){
Realpart = rnorm(1,0, sqrt(V[k]/2)); 
Imagpart = rnorm(1,0, sqrt(V[k]/2));
PosCoeff[k] = complex(real=Realpart, imag=Imagpart);
}

RandCoeff[(T/2):T] = PosCoeff;
RandCoeff[1:(T/2-1)] = Conj(PosCoeff[(T/2):-1:2]);


### plot the real part of the coefficents
par(mfrow=c(1,2)); 
plot(Re(RandCoeff), type="l");
plot(Im(RandCoeff), type="l");

### plot the squared coefficients
par(mfrow=c(1,1)); 
Pgram = (abs(RandCoeff))^2;
plot(Pgram, type="l")

### Compare this with the variance function 
par(mfrow=c(1,1)); 
ymax = max(Pgram, V); ymin = min(Pgram, V); 
plot((1:(T/2+1)), V[1:(T/2+1)], ylim=c(ymin, ymax), xlab="Frequency", ylab="Var", col=1); 
lines(Pgram[(T/2):T], col=2)

### The time series 
tsdata = FB%*%RandCoeff;
plot(Im(tsdata))
tsdata = Re(tsdata);
par(mfrow=c(1,1)); 
plot(tsdata, type="l")

############################################################################
############################################################################
############### Estimating the variance decomposition ######################
############### from a given time series ###################################

########## Example 1
T=1024;
ts1 = arima.sim(T, model=list(ar=c(0.9)));
par(mfrow=c(1,1)); plot(ts1, type="l");

## Compute the Fourier coefficients 
fft1 = (1/sqrt(T))*fft(ts1); 

## Squared Coefficients
sqcoef = (abs(fft1))^2;

## compare energy in time and frequency domains 
sum(ts1^2)
sum(sqcoef)

## Periodograms
Pgram1 = (1/(2*pi))*sqcoef;
par(mfrow=c(1,1)); plot(Pgram1, type="l")

Pgram1x = Pgram1[1:(T/2+1)];
par(mfrow=c(1,1)); plot(Pgram1x, type="l")

########### Example 2

T=1024;
ts1 = arima.sim(T, model=list(ar=c(-0.9)));
par(mfrow=c(1,1)); plot(ts1, type="l");

## Compute the Fourier coefficients 
fft1 = (1/sqrt(T))*fft(ts1); 

## Squared Coefficients
sqcoef = (abs(fft1))^2;

## compare energy in time and frequency domains 
sum(ts1^2)
sum(sqcoef)

## Periodograms
Pgram2 = (1/(2*pi))*sqcoef;
par(mfrow=c(1,1)); plot(Pgram2, type="l")

Pgram2x = Pgram2[1:(T/2+1)];
par(mfrow=c(1,1)); plot(Pgram2x, type="l")

############ Example 3

T=1024;
ts1 = rnorm(T,0,1);
par(mfrow=c(1,1)); plot(ts1, type="l");

## Compute the Fourier coefficients 
fft1 = (1/sqrt(T))*fft(ts1); 

## Squared Coefficients
sqcoef = (abs(fft1))^2;

## compare energy in time and frequency domains 
sum(ts1^2)
sum(sqcoef)

## Periodograms
Pgram3 = (1/(2*pi))*sqcoef;
par(mfrow=c(1,1)); plot(Pgram3, type="l")

Pgram3x = Pgram3[1:(T/2+1)];
par(mfrow=c(1,1)); plot(Pgram3x, type="l")


############################################################################
############################################################################
############### Spectrum of stationary processes ###########################
############################################################################

########## Example 1
T = 1024;
spec1 = armaspec(T, arcoeff=c(0.9), macoeff=c(0), vare=1);
plot(spec1, type="l")

ymax = max(spec1, Pgram1x); ymin=min(spec1, Pgram1x);
plot(Pgram1x, type="l", col=1, ylim=c(ymin, ymax)); 
lines(spec1, col=2, lwd=2);

########## Example 2
T = 1024;
spec2 = armaspec(T, arcoeff=c(-0.9), macoeff=c(0), vare=1);
plot(spec2, type="l")

ymax = max(spec2, Pgram2x); ymin=min(spec2, Pgram2x);
plot(Pgram2x, type="l", col=1, ylim=c(ymin, ymax)); 
lines(spec2, col=2, lwd=2);

########## Example 3
T = 1024;
spec3 = armaspec(T, arcoeff=c(0), macoeff=c(0), vare=1);
plot(spec3, type="l")

ymax = max(spec3, Pgram3x); ymin=min(spec3, Pgram3x);
plot(Pgram3x, type="l", col=1, ylim=c(ymin, ymax)); 
lines(spec3, col=2, lwd=2);


################################################################################
################################################################################
############### Estimating by kernel smoothing #################################
################################################################################

#### Example 1 

smper = periodsmoother(Pgram1, 1);
plot((1:(T/2+1)), Pgram1x, type="l", xlab="Freq", ylab="");
lines(smper, col=2)

ymax = max(spec1, smper); ymin=min(smper, spec1);
plot((2:(T/2+1)), smper, type="l", ylim=c(ymin, ymax), xlab="Freq", ylab="");
lines(spec1, col=2)

#################################################################################
################### Filtering ###################################################
#################################################################################

T=1024;
ts1 = arima.sim(T, model=list(ar=c(0.9)));
ts2 = arima.sim(T, model=list(ar=c(-0.8)));

par(mfrow=c(1,1)); plot(ts1, type="l");
par(mfrow=c(1,1)); plot(ts2, type="l");

Y = ts1 + ts2;
par(mfrow=c(1,1)); plot(Y, type="l");

fft1 = (1/sqrt(T))*fft(Y);
periodY = (1/2*pi)*((abs(fft1))^2);
periodYx = Re(periodY[1:(T/2)]);
plot(periodYx, type="l");


## Suppose want to remove the high frequencies
  
powertfout = powertf(T, c(1/5,1/5), c(1/5), c(1/5,1/5));
plot(powertfout, type="l")

filtspecest = powertfout*periodYx;
plot(filtspecest[1:(T/2)], type="l", col=1);
lines(periodYx[1:(T/2)], col=2)