Consider the vector AR(3) process (“~ 20 Hz β-network”)
$$ \begin{align} X_{1,t} &= 0.95\sqrt(2)X_{1,t−1} − 0.9025X_{1,t−2} + ε_{1,t} \\ X_{2,t} &= 0.5X_{1,t−2} + ε_{2,t} \\ X_{3,t} &= −0.4X_{1,t−3} + ε_{3,t} , X_{4,t} = −0.5X_{1,t−2} + 0.25\sqrt(2)X_{4,t−1} + 0.25\sqrt(2)X_{5,t−1} + ε_{4,t} \\ X_{5,t} &=−0.25\sqrt(2)X_{4,t−1}+0.25\sqrt(2)X_{5,t−1}+ε_{5,t} \end{align} $$where $X_{t}$ is a 5-dimensional column vector. The noise vector is Gaussian white with diagonal covariance matrix, where the diagonal entries correspond to [0.6, 0.5, 0.3, 0.3, 0.6]. (Note: a diagonal matrix is used here for simplicity. We could explore instantaneous noise correlations with a nondiagonal matrix.)
(1 point) Using the script provided in the warm-up code, sample 200 trials, each 3-second long with a sampling rate of 200 Hz. The initial condition can be zero mean unit variance Gaussian. Remember to throw way the first 200 samples or so to remove the transient. Plot a single trial showing all of the 5 channels. For sampling, we define the coefficient matrix in Matlab as (M = 5 is the dimension of the AR process of order = 3)
p = 3 #order of AR(p)
dim = 5 #number of random variables
A = zeros(Float64, dim, dim, p)
a = sqrt(2)
A[1,1,1] = 0.95*a
A[1,1,2] = -0.9025
A[2,1,2] = 0.5
A[3,1,3] = -0.4
A[4,1,2] = -0.5
A[4,4,1] = 0.25*a
A[4,5,1] = 0.25*a
A[5,4,1] = -0.25*a
A[5,5,1] = 0.25*aAdapt the script provided in the warm-up code to estimate the spectral matrix based on the multitaper method. Use the provided function “multitaperSpectrum.m” and a 5 Hz bandwidth.
(3 points) Compute and plot the power spectrum, spectral coherence and partial spectral coherence based on the estimated spectral matrix. (Your results should replicate the figure shown in the Lecture 5.)
In [1]:
using Distributions
using PlotlyJS
In [2]:
# AR matrices
p = 3 #order of AR(p)
dim = 5 #number of random variables
T=3# Trial duration in seconds
ntrials=200
Fs=200# sampling rate in per second
dt=1/Fs
burnin=1000 # for transient removal
df = 1/T; #frequency resolution
N=T*Fs
Out[2]:
In [3]:
# Use companion matrix to check for stability
A = zeros(Float64, dim, dim, p)
a = sqrt(2)
A[1,1,1] = 0.95*a
A[1,1,2] = -0.9025
A[2,1,2] = 0.5
A[3,1,3] = -0.4
A[4,1,2] = -0.5
A[4,4,1] = 0.25*a
A[4,5,1] = 0.25*a
A[5,4,1] = -0.25*a
A[5,5,1] = 0.25*a
#Companion matrix
m,n,p = size(A)
pn = (p-1)*m
Ac = [reshape(A,m,p*n); eye(pn) zeros(pn,m)];
In [4]:
#stationary?
D, V = eig(Ac)
#spectracl radius
ρ = maximum(abs.(D))
if ρ < 1.0
println("Process is Stationary")
else
println("Process is NOT Stationary")
end
println("Spectral Radius ρ(Ac) = ", ρ)
In [5]:
# Sampling from the VAR(3) model
# The noise covariance (Using DIAGONAL COVARIANCE MATRIX FOR SIMPLICITY ...)
ntrials = 200
Σ_diag = [0.6, 0.5, 0.3, 0.3, 0.6]
μ = zeros(Float64, dim)
d = MvNormal(μ, Σ_diag)
Out[5]:
In [6]:
X = zeros(Float64, dim, N+burnin, ntrials)
for r=1:ntrials
#Multivariate GWN sequence of dimension (dim x cases)
E=rand(d, N+burnin)
#Intial state - first "p" times
X[:,1:p,r]=E[:,1:p]
#Iterate through time
for k=p+1:N+burnin
x_k=zeros(Float64, dim)
#autoregression
for j=1:p
x_k=x_k+A[:,:,j]*X[:,k-j,r] #5-dim
end
X[:,k,r]=x_k + E[:,k]
end
end
#remove burn-in
X=X[:,burnin+1:end,:];
In [7]:
#Plot the path
plots = Vector{typeof(plot())}(dim)
for d=1:dim
plots[d] = plot(0:dt:dt*N, X[d,:,1], line_width=0.8, name="X$d");
end
all_plots = [plots...]
all_plots.plot.layout["title"] = "VAR(3)"
all_plots.plot.layout["xaxis5_title"] = "time"
all_plots
Out[7]:
In [8]:
using DSP
using FFTW
include("transform_stats.jl")
include("multitaper.jl")
Out[8]:
In [9]:
#need to permute dims because multitaper expects matrix of dims = samples x channels x trials
A = permutedims(X, [2,1,3])
fNyquist=Fs/2;
freqrange = 0:df:fNyquist#frequency grid
bandWidth = 5
NFFT=N
(xs, s) = multitaper(A, [CrossSpectrum(), PowerSpectrum()], Fs, tapers = dpss(size(A, 1), bandWidth), nfft = N);
In [10]:
# Coompute the coherence and partial coherence
pairs = allpairs(dim)
cop = zeros(Float64, length(freqrange), size(pairs, 2))
co = zeros(Float64, length(freqrange), size(pairs, 2))
for i = 1:size(pairs, 2)
ch1 = pairs[1, i]
ch2 = pairs[2, i]
i1=setdiff(1:dim,ch1)
i2=setdiff(1:dim,ch2)
for f=1:length(freqrange)
S = eye(Complex{Float64}, dim).*s[f,:]
for j = 1:size(pairs, 2)
S[pairs[1, j], pairs[2, j]] = xs[f, j]
S[pairs[2, j], pairs[1, j]] = conj(xs[f, j])
end
M12 = det(S[i1,i2])
M11 = det(S[i1,i1])
M22 = det(S[i2,i2])
co[f, i] = abs(xs[f, i])^2/(s[f, ch1]*s[f, ch2])
cop[f, i] = abs(M12)/sqrt(abs(M11)* abs(M22))
end
end
In [11]:
#---------------- Plot Coherence and Partial Coherence---------------------
nplots_c = size(pairs, 2)
plots_c = Vector{typeof(plot())}(nplots_c)
for i = 1:nplots_c
ch1 = pairs[1, i]
ch2 = pairs[2, i]
c_trace = scatter(x = collect(freqrange), y=co[:, i], line_width=0.8, name="X$i", marker_color="red")
cp_trace = scatter(x=collect(freqrange), y=cop[:, i], line_width=0.8, name="X$i", marker_color="green");
plots_c[i] = plot([c_trace, cp_trace]);
end
#---------------- Plot Power Spectrum---------------------
plots_s = Vector{typeof(plot())}(dim)
for i=1:dim
plots_s[i] = plot(collect(freqrange), s[:, i], line_width=0.8, name="X$i", marker_color="blue");
end
all_plots = [plots_s[1] plots_c[1] plots_c[2] plots_c[3] plots_c[4];
plots_c[1] plots_s[2] plots_c[5] plots_c[6] plots_c[7];
plots_c[2] plots_c[5] plots_s[3] plots_c[8] plots_c[9];
plots_c[3] plots_c[6] plots_c[8] plots_s[4] plots_c[10];
plots_c[4] plots_c[7] plots_c[9] plots_c[10] plots_s[5]]
all_plots.plot.layout["title"] = "Power(blue) , Square-Coherence(red), Partial-Coherence (green)"
all_plots.plot.layout["showlegend"] = false
all_plots
Out[11]:
Based on the sample from the model generated in exercise (1),
In [12]:
include("var_maxent.jl");
In [13]:
# Estimating VAR(p) from data ===================
df_est=0.001#frequency resolution to be used in the parametric AR estimation
nF=Int(round(1/df_est+1))#number of frequencies evaluated in the parametric AR estimation
Ah,SIGMAh,Eh = var_maxent(X, p)
DSIG = det(SIGMAh)# residuals covariance matrix determinant
if DSIG <= 0
warn("Residuals covariance not positive definite")
end
M=N*ntrials
L= -(M/2)*log(DSIG) #max loglikelihood
aic = -2*L + 2*p*dim^2*(M/(M-p-1)) # Note AIC without correction = -2*L + 2*order*xDim^2
bic = -2*L + p*log(M)
# Transfer function H(f) and spectral matrix S(f) from fitted model
# F=collect(0:nF-1)*df*Fs/2
delta=1
F=linspace(0,0.5,nF);
j=0
Sh = zeros(Complex{Float64}, dim, dim, nF)
for f in F
j=j+1
H = eye(dim) # identity matrix
for m=1:p
H=H-Ah[:,:,m]*exp(-1im*m*2*pi*f)
end
H = inv(H)
#S(:,:,j) = H*SIGMAh*ctranspose(H)
Sh[:,:,j] = H*SIGMAh*H'
#Note: MATLAB (julia - should also do the same through multiple dispatch)
#1. the transpose in will also recognize the complex data and perform the conjugate transpose
#2. however, that .' implements only the transpose, not the conjugate transpose
end
Sh = 2 * delta^2 * 1/(N*delta) * Sh #Normalize the spectral matrix
for k = 1:dim
Sh[k,k,:] = real(Sh[k,k,:])
end #make sure the diagonal has real numbers (power spectrum)
# F=F*Fs# Change frequency to Hz;
In [14]:
# To compute the partial coherence =======================
# S is the spectral matrix S(channel j, channel k, frequency), j,k = 1, 2, ... xDim
Ch = zeros(Float64, nF, size(pairs, 2))
CPh = zeros(Float64, nF, size(pairs, 2))
Ph = zeros(Float64, nF, dim)
for i = 1:size(pairs, 2)
j = pairs[1, i]
k = pairs[2, i]
i1=setdiff(1:dim,j)
i2=setdiff(1:dim,k)
for f=1:nF
Mjk = det(Sh[i1,i2,f])
Mjj = det(Sh[i1,i1,f])
Mkk = det(Sh[i2,i2,f])
Ch[f, i] = abs(Sh[j,k,f])^2/real(Sh[j,j,f]* Sh[k,k,f])
CPh[f, i] = abs(Mjk)/real(sqrt(Mjj*Mkk))#to avoid numerical issues, force it to be real
end
end
for j=1:dim
for f=1:nF
Ph[f,j] = abs(Sh[j,j,f])
end
end
In [15]:
#---------------- Plot Coherence and Partial Coherence---------------------
nplots_c = size(pairs, 2)
plots_c = Vector{typeof(plot())}(nplots_c)
for i = 1:nplots_c
ch1 = pairs[1, i]
ch2 = pairs[2, i]
c_trace = scatter(y=Ch[:, i], line_width=0.8, name="X$i", marker_color="red")
cp_trace = scatter(y=CPh[:, i], line_width=0.8, name="X$i", marker_color="green");
plots_c[i] = plot([c_trace, cp_trace]);
end
#---------------- Plot Power Spectrum---------------------
plots_s = Vector{typeof(plot())}(dim)
for i=1:dim
plots_s[i] = plot(Ph[:, i], line_width=0.8, name="X$i", marker_color="blue");
end
all_plots = [plots_s[1] plots_c[1] plots_c[2] plots_c[3] plots_c[4];
plots_c[1] plots_s[2] plots_c[5] plots_c[6] plots_c[7];
plots_c[2] plots_c[5] plots_s[3] plots_c[8] plots_c[9];
plots_c[3] plots_c[6] plots_c[8] plots_s[4] plots_c[10];
plots_c[4] plots_c[7] plots_c[9] plots_c[10] plots_s[5]]
all_plots.plot.layout["title"] = "Power(blue) , Square-Coherence(red), Partial-Coherence (green)"
all_plots.plot.layout["showlegend"] = false
all_plots
Out[15]:
In [ ]: