Frequency Domain Connectedness

Conenctendess

VAR:

\[ \begin{equation}\mathbf{x}_t = \mathbf{\Phi}_1 \mathbf{x}_{t-1} + \mathbf{\Phi}_2 \mathbf{x}_{t-2} + \cdots + \mathbf{\Phi}_p \mathbf{x}_{t-p} + \boldsymbol{\epsilon}_t,\end{equation} \]

with \(\Phi_1,...,\Phi_p\) coefficient matrices, and \(\epsilon_t\) being white noise with (possibly non-diagonal) covariance matrix \(\Sigma\).

Expression in Lag Operator:

\[ \begin{equation}\mathbf{\Phi}(L)\mathbf{x}_t = \boldsymbol{\epsilon}_t\end{equation} \]

VMA Representation (MA(\(\infty\))):

\[ \begin{equation}\mathbf{x}_t = \Psi(L) \boldsymbol{\epsilon}_t,\end{equation} \]

where \(\Phi(L)\) matrix of infinite lag polynomials.

Joint Connectedness

\[ \begin{align*} \mathcal{S}^{jnt}_{\bullet \to i} = \frac{E\left(\xi_{i,t}^2(H)\right) - E\left[\xi_{i,t}(H) - E(\xi_{i,t}(H) \mid \epsilon_{\forall j \neq i, t+1, \ldots, \epsilon_{\forall j \neq i, t+H}} )\right]^2} {E\left(\xi_{i,t}^2(H)\right)} \\ = \frac{ \sum_{h=0}^{H-1} e_i' A_h \Sigma_\epsilon M_i (M_i' \Sigma_\epsilon M_i)^{-1} M_i' \Sigma_\epsilon A_h' e_i }{ \sum_{h=0}^{H-1} e_i' A_h \Sigma_\epsilon A_h' e_i } \end{align*} \]where \(A_h\) are the VMA coefficient matrices,
\(\Sigma\) the covariance matrix of white noise shocks,
\(M_i\) the elimination matrix, a rectangular matrix that equals \(I\) with the \(i\)th column eliminated
\(e_i\) is a selection vector containing all zeroes except for a 1 in the \(i\)th element

Calculation in R

A = Wold(Phi[,,ij], nfore)  # the VMA coefficient matrices
Xi_h = array(0,dim=c(k,k,nfore))
for (h in 1:nfore) {
  Xi_h[,,h] = A[,,h]%*%Sigma[,,ij]%*%t(A[,,h]) # calculated Xi at each h
}
Xi = rowSums(Xi_h, dims=2) # sum them along THIRD dimension to form Xi  (note: because this is a row sum, dims=2, actually sums along the third dimension)
I_K = diag(1,nrow=k,ncol=k)

#Calculate the elimination matrices.
#Usually denoted as a KxK-1 matrix M_i, here it is an array where M[,,1]=M_1, and in general M[,,i]=M_i.
M = array(0,dim=c(k,k-1,k))
for (i in 1:k){
  M[,,i] = I_K[,-i]  #calculate the elimination matrices
}

#Calculate the joint total spillovers from all others to variable i (S_jnt)
#calculate the numerator of S_jnt
S_jnt_numerator_h = array(0,dim=c(k,nfore))
for (i in 1:k){
  for (h in 1:nfore){
    S_jnt_numerator_h[i,h] = I_K[i,]%*%A[,,h]%*%Sigma[,,ij]%*%M[,,i]%*%(solve(t(M[,,i])%*%Sigma[,,ij]%*%M[,,i]))%*%t(M[,,i])%*%Sigma[,,ij]%*%t(A[,,h])%*%I_K[,i]     #calculate the numerator of S_jnt at each h
  }
}
S_jnt_from = array(0,dim=c(k))
for (i in 1:k){
  S_jnt_from[i] = 100*sum(S_jnt_numerator_h[i,])/Xi[i,i]
}

Frequency Connectedneess

Frequency Response Function: 1

  • the Fourier transform of the impulse response \(\Phi_h\)

\[ \begin{equation}\Psi\left(e^{-i\omega}\right) = \sum_{h} e^{-i\omega h} \Psi_h \tag*{[BK18]}\end{equation} \]

Special Density Matrix (Power Spectral Density):

\[ \begin{equation}S_{\mathbf{x}}(\omega) = \sum_{h=-\infty}^{\infty} E\left(\mathbf{x}_t \mathbf{x}_{t-h}'\right) e^{-i \omega h}= \Psi\left(e^{-i\omega}\right)\Sigma\Psi'\left(e^{+i\omega}\right).\tag*{[BK18]}\end{equation} \]

The Generalized Causation Spectrum

\[ \begin{equation} \left(\tilde{\mathbf{f}}(\omega)\right)_{j,k} \equiv \frac{ \sigma_{kk}^{-1} \left| \left( \Psi(e^{-i\omega}) \Sigma \right)_{j,k} \right|^2 }{ \left( \Psi(e^{-i\omega}) \Sigma \Psi'(e^{+i\omega}) \right)_{j,j} }, \end{equation} \] The scaled generalized variance decomposition

  • on the frequency band \(d = (a, b) : a, b \in (-\pi, \pi),\, a < b\):

\[ \begin{equation} (\tilde{\theta}_d)_{j,k} = (\theta_d)_{j,k} \Big/ \sum_{k} (\theta_{\infty})_{j,k}, \end{equation} \]

In sample…

\[ \begin{equation} \left( \hat{\theta}_d \right)_{j,k} = \sum_{\omega} \hat{\Gamma}_j(\omega) \left( \left( \hat{\tilde{f}}(\omega) \right)_{j,k} \right), \end{equation} \]

where

\[ \begin{equation} \left( \hat{\tilde{f}}(\omega) \right)_{j,k} \equiv \frac{ \hat{\sigma}_{kk}^{-1} \left| \left( \hat{\Psi}(\omega) \hat{\Sigma} \right)_{j,k} \right|^2 }{ \left( \hat{\Psi}(\omega) \hat{\Sigma} \hat{\Psi}'(\omega) \right)_{j,j} } \end{equation} \]

is the estimated generalized causation spectrum, and

\[ \begin{equation} \hat{\Gamma}_j(\omega) = \frac{ \left( \hat{\Psi}(\omega) \hat{\Sigma} \hat{\Psi}'(\omega) \right)_{j,j} }{ \left( \Omega \right)_{j,j} } \end{equation} \]

is the estimate of the weighting function, where \(\Omega = \sum_{\omega} \hat{\Psi}(\omega) \hat{\Sigma} \hat{\Psi}'(\omega).\)

R Implementation

1fftir = lapply(irf, function(i) apply(i, 2, fft))
fftir = lapply(1:(nfore + 1), function(j) sapply(fftir, function(i) i[j, ]))
2denom = diag(Re(Reduce("+", lapply(fftir, function(i) i %*% Sigma %*% t(Conj(i))/nfore)[range])))
if (generalized) {
  irf_ = lapply(fftir, function(i) t(abs(i %*% Sigma) %*% solve(diag(sqrt(diag(Sigma))))))
  IRF = array(unlist(irf_), c(k,k,nfore+1))
3  irf_ = lapply(fftir, function(i) (abs(i %*% Sigma)))
  IRF_ = array(unlist(irf_), c(k,k,nfore+1))
4  enum_ = IRF_^2/(nfore + 1)
  tab = list()
  for (i in 1:dim(enum_)[3]) {
    fevd_ = enum_[,,1]
    for (j in 1:dim(enum_)[2]) {
5      fevd_[j,] = enum_[j,,i] / (denom[j] * diag(Sigma))
    }
    tab[[i]] = t(fevd_)
  }
  tot = apply(Reduce("+", tab[range]), 2, sum)
6  FEVD = lapply(tab, function(i) t(i)/tot)
}
1
perform FFT on IRF
2
the denominator of \(\hat{\Gamma}_j(\omega)\) (* the total version, rather than the WTH version)
3
the absolute value part of the numerator
4
squaring the absolute value part of the numerator
5
diag(Sigma) -> the \(sigma^-1_kk\) part of the numerator
6
the scaled generalized variance decomposition

The weighting Function

  • to convert from within-frequency causation to a natural decomposition of variance decompositions to frequencies

\[ \begin{equation} \Gamma_j(\omega) = \frac{ \left(\Psi(e^{-i\omega}) \Sigma \Psi'(e^{+i\omega})\right)_{j,j} }{ \frac{1}{2\pi} \int_{-\pi}^{\pi} \left(\Psi(e^{-i\lambda}) \Sigma \Psi'(e^{+i\lambda})\right)_{j,j} \, d\lambda }, \end{equation} \]

Relationship with Time Domain Connectedness

\[ \begin{equation} (\theta_{\infty})_{j,k} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \Gamma_j(\omega)\, (\tilde{\mathbf{f}}(\omega))_{j,k} \, \mathrm{d}\omega. \end{equation} \]

How does the freuqency partitioning works?

The freuqency partitioning is specified as: (Example)

data("cgg2022")
partition = c(pi+0.00001, pi/5, 0)
dca = ConnectednessApproach(cgg2022, 
                            model="TVP-VAR",
                            connectedness="Frequency",
                            nlag=1,
                            nfore=100,
                            window.size=200,
                            VAR_config=list(TVPVAR=list(kappa1=0.99, kappa2=0.99, prior="BayesPrior")),
                            Connectedness_config = list(
                              FrequencyConnectedness=list(partition=partition, generalized=TRUE, scenario="ABS")
                            ))

In FrequencyConnectedness.R at ConnectednessApproach, the specified frequency partitions are then converted into a group of indexes. For a simple case of 2 partitions of c(pi+0.0001, pi/5, 0)

Result is

[[1]]
 [1] 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
[46] 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92

[[2]]
 [1]   1   2   3   4   5   6   7   8   9  10  93  94  95  96  97  98  99 100 101
  new_p = frequencyConnectedness::getPartition(partition, nfore)
  range = sort(unique(do.call(c, new_p)))
  
  date = as.character(date)
  TCI = array(0, c(t,interval), dimnames=list(date, period_names))
  PCI = INFLUENCE = CT = NPDC = array(0, c(k, k, t, interval), dimnames=list(NAMES, NAMES, date, period_names))
  NET = FROM = TO = array(0, c(t, k, interval), dimnames=list(date, NAMES, period_names))
  NPT = array(0, c(t, k, interval), dimnames=list(date, NAMES, period_names))
  PCI = INFLUENCE = array(0, c(k, k, t, interval), dimnames=list(NAMES, NAMES, date, period_names))
  pb = progress_bar$new(total=t)
  for (i in 1:t) {
    decomp = FEVD(Phi=Phi[,,i], Sigma=Sigma[,,i], nfore=nfore, generalized=generalized, type="frequency", range=range)$FEVD
    for (ij in 1:length(decomp)) {
      rownames(decomp[[ij]]) = colnames(decomp[[ij]]) = 1:ncol(Sigma)
    }
    tables = lapply(new_p, function(j) Reduce('+', decomp[j]))
    for (j in 2:interval) {
      if (scenario=="ABS") {
        dca = ConnectednessTable(tables[[j-1]])
        CT[,,i,j] = dca$FEVD
        TO[i,,j] = dca$TO
        FROM[i,,j] = dca$FROM
        NET[i,,j] = dca$NET
        NPDC[,,i,j] = dca$NPDC
        INFLUENCE[,,i,j] = dca$INFLUENCE
        NPT[i,,j] = dca$NPT
        if (corrected) {
          TCI[i,j] = dca$cTCI
        } else {
          TCI[i,j] = dca$TCI
        }
      } else if (scenario=="WTH") {
        dca = ConnectednessTable(tables[[j-1]]/sum(sum(tables[[j-1]]))*k)
        CT[,,i,j] = dca$FEVD
        TO[i,,j] = dca$TO
        FROM[i,,j] = dca$FROM
        NET[i,,j] = dca$NET
        NPDC[,,i,j] = dca$NPDC
        INFLUENCE[,,i,j] = dca$INFLUENCE
        NPT[i,,j] = dca$NPT
        if (corrected) {
          TCI[i,j] = dca$cTCI
        } else {
          TCI[i,j] = dca$TCI
        }
      }
    }
getIndeces <- function(n.ahead, up, down) {
    space <- (0:floor(n.ahead/2))/((floor(n.ahead/2)))*pi

    lb <- space >= down
    ub <- space < up # That's why the "up" has to be specified with like pi+0.00001 ....
    output = (lb & ub)*1
    if (n.ahead%%2 == 0) {
        output <- c(output, rev(output[2:(length(output)-1)]))
    }
    else {
        output <- c(output, rev(output[2:length(output)]))
    }
    return(which(output==1))
}

getPartition <- function(partition, n.ahead) {
    if (!all(sort(partition, decreasing = T)==partition)) {
        stop("The bounds must be in decreasing order.")
    }
    part <- lapply(1:(length(partition)-1), function(i) getIndeces(n.ahead+1, partition[i], partition[i+1]))
    if (length(unique(do.call(c, part))) != (n.ahead + 1)) {
            warning("The selected partition does not cover the whole range.")
    }
    if (any(sapply(part, length)==0)) {
        sprintf("The n.ahead steps does not allow to infer anything about the following interval (%f, %f).", partition[which(sapply(part, length)==0)], partition[which(sapply(part, length)==0)+1])
        stop("Change the partition.")
    }
    return(part)
}

An example

  • Property of (discrete REAL) FFT
r_n <- 10
r <- rnorm(r_n)
fftr <- fft(r)

matrix(fftr, r_n, 1)
                        [,1]
 [1,]  -3.2571483+0.0000000i
 [2,]  -0.6999804-0.5103578i
 [3,]   1.8202231-2.5015829i
 [4,]  -2.3051612+1.1744117i
 [5,]  -1.0892069+2.0296636i
 [6,] -10.3208437+0.0000000i
 [7,]  -1.0892069-2.0296636i
 [8,]  -2.3051612-1.1744117i
 [9,]   1.8202231+2.5015829i
[10,]  -0.6999804+0.5103578i
exp(-2 * pi * 1i * 0.5)
[1] -1-0i
# fc is a function to evaluate the *frequency component* of value or series `v` at specified frequency `f`
fc <- function(v, f) {
  t <- 1:length(v)
  sum(v * exp(-2 * pi * 1i * f * (t - 1)))
}

# Check fft formula
?fft
starting httpd help server ... done
fc(0.25, 1/2)
[1] 0.25+0i
fc(r, 1/2)
[1] -10.32084-0i
# a *Fourier Transform* return the frequency component at different frequency from 0 to (n - 1) / n
ft <- function(v) sapply((0:(length(v)-1))/length(v), function(f) fc(v, f))

ft(r)
 [1]  -3.2571483+0.0000000i  -0.6999804-0.5103578i   1.8202231-2.5015829i
 [4]  -2.3051612+1.1744117i  -1.0892069+2.0296636i -10.3208437-0.0000000i
 [7]  -1.0892069-2.0296636i  -2.3051612-1.1744117i   1.8202231+2.5015829i
[10]  -0.6999804+0.5103578i
fft(r)
 [1]  -3.2571483+0.0000000i  -0.6999804-0.5103578i   1.8202231-2.5015829i
 [4]  -2.3051612+1.1744117i  -1.0892069+2.0296636i -10.3208437+0.0000000i
 [7]  -1.0892069-2.0296636i  -2.3051612-1.1744117i   1.8202231+2.5015829i
[10]  -0.6999804+0.5103578i
# an illustration of the unit circle around exp(-2 * pi * 1i * x)
plot(sapply((0:99)/100, function(x) exp(- 2 * pi * 1i * x)), asp = 1)

Hermetian Symmetry

  • Note the important property of the “Hermetian symmetry” - the second half of the fft series on real number are the complex conjugate of the first half
matrix(fftr, length(fftr), 1)
                        [,1]
 [1,]  -3.2571483+0.0000000i
 [2,]  -0.6999804-0.5103578i
 [3,]   1.8202231-2.5015829i
 [4,]  -2.3051612+1.1744117i
 [5,]  -1.0892069+2.0296636i
 [6,] -10.3208437+0.0000000i
 [7,]  -1.0892069-2.0296636i
 [8,]  -2.3051612-1.1744117i
 [9,]   1.8202231+2.5015829i
[10,]  -0.6999804+0.5103578i

Why the “Hermetian Symmetry”?

The “basis function” \(e^{-2\pi i n k/N}\) or \(e^{-2\pi i (t - 1) f},\) where \(f = k/T, k = 0/T, 1/T, ... , (T-1)/T\) (as expressed in the R way) in the fft formula: \(\sum^T_{t=1}z[k] \cdot e^{-2\pi i (t - 1)f}\)

if \(1/2 < f < 1\) (per sampling rate, e.g. day), the Nyquist Frequency,

the number \(f\) can be expressed \(1 - p\), where \(p\) is a positive number with \(p < 1/2\).

Then the basis function can be written as \(e^{-2\pi i (t - 1) f} = e^{-2\pi i (t - 1) (1 - p)}\). Due to the periodicity of complex exponentials \(e^{-2\pi i x} = 1, \forall x \in \mathbb{Z}\), \(e^{-2\pi i (t - 1) (1 - p)} = 1 \cdot e^{-2\pi i (t-1) (-p)} = e^{2\pi i (t-1) p}\)

We have

\[ e^{-2\pi i (t-1) f} = e^{2\pi i (t-1) p} \]

That is, if \(1/2 , f < 1\), \(e^{-2\pi i (t-1) f}\) is just the complex conjugate of \(e^{2\pi i (t-1) p}, p = 1 - f\)

More on complex conjugate:

With Euler’s formula \(e^{i\theta} = \cos \theta + i \sin \theta\),

The switch of side in the exponent yield \(e^{-i\theta} = -\cos \theta - i \sin \theta = cos \theta - i \sin \theta\), which is the complex conjugate of e^{i}

Negative Frequency

Due to the Hermetian Symmetry (and the fact that \(e^{i\theta} = Conj(e^{-i\theta})\)), the output of FFT on real numbers over the frequency \((0, 2\pi)\) or \((-\pi, \pi)\) are just a shift of position. The sum of results over these 2 ranges / frequency bands are the same (BK18 uses the \((-\pi, \pi) convention\)). The result of the fft of a negative frequency \(f\) is the Complex conjugate of the positive frequency counterpart \(|f|\)

exp(-2 * pi * 1i * (0:10)/11)
 [1]  1.0000000+0.0000000i  0.8412535-0.5406408i  0.4154150-0.9096320i
 [4] -0.1423148-0.9898214i -0.6548607-0.7557496i -0.9594930-0.2817326i
 [7] -0.9594930+0.2817326i -0.6548607+0.7557496i -0.1423148+0.9898214i
[10]  0.4154150+0.9096320i  0.8412535+0.5406408i
exp(-2 * pi * 1i * (-5:5)/11)
 [1] -0.9594930+0.2817326i -0.6548607+0.7557496i -0.1423148+0.9898214i
 [4]  0.4154150+0.9096320i  0.8412535+0.5406408i  1.0000000+0.0000000i
 [7]  0.8412535-0.5406408i  0.4154150-0.9096320i -0.1423148-0.9898214i
[10] -0.6548607-0.7557496i -0.9594930-0.2817326i
sum(exp(-2 * pi * 1i * (0:10)/11)) # they are the same, after correcting computational errors
[1] -3.885781e-16-3.885781e-16i
sum(exp(-2 * pi * 1i * (-5:5)/11))
[1] 5.551115e-16+0i

What happen in the Estimation of Frequency Connectedness, step by step

pacman::p_load(vars, ConnectednessApproach, xts)

n = 1000
nfore = 100
d <- xts(cbind(rnorm(n),rnorm(n)), order.by = seq.Date(from = as.Date("2023-01-01"), by = "day", length.out = n))
var_model <- VAR(d, configuration = list(nlag=1))

irf_ = IRF(var_model$B, var_model$Q, nfore=nfore, orth=FALSE)$irf

Sigma = var_model$Q[,,1] # # hardcoded t at Sigma
k = ncol(Sigma)

partition = c(pi+0.00001, pi/5, 0)
new_p = frequencyConnectedness::getPartition(partition, nfore)
range = sort(unique(do.call(c, new_p)))

fftir = lapply(irf_, function(i) apply(i, 2, fft))
fftir = lapply(1:(nfore + 1), function(j) sapply(fftir, function(i) i[j, ]))

# irf__ = lapply(1:(nfore + 1), function(j) sapply(irf_, function(i) i[j, ]))
# Reduce("+", lapply(irf__, function(i) i %*% Sigma %*% t(Conj(i))/nfore))

denom = diag(Re(Reduce("+", lapply(fftir, function(i) i %*% Sigma %*% t(Conj(i))/nfore)[range])))

irf_ = lapply(fftir, function(i) (abs(i %*% Sigma)))
IRF_ = array(unlist(irf_), c(k,k,nfore+1))
enum_ = IRF_^2/(nfore + 1)
tab = list()
for (i in 1:dim(enum_)[3]) {
  fevd_ = enum_[,,1]
  for (j in 1:dim(enum_)[2]) {
    fevd_[j,] = enum_[j,,i] / (denom[j] * diag(Sigma))
    
    print(fevd_[j,])
  }
  tab[[i]] = t(fevd_)
}
tot = apply(Reduce("+", tab[range]), 2, sum)
FEVD_ = lapply(tab, function(i) t(i)/tot)

tables = lapply(new_p, function(j) Reduce('+', FEVD_[j]))

Why is the interpretation of “1-5 days” for the frequency band \((pi, pi/5)\) probably wrong?

  • The Nyquist frequency of a discrete time series \(f_n\) is \(1/2\) (half of the sampling frequency), which is the highest frequency reliable measurable in the data. The events at frequency \(f > 1/2\), i.e. < 2 days a cycle, therefore couldn’t be measurable correctly.
  • \((pi, pi/5)\) does not correspond to a frequency band of 5 days a cycle or less, using the getIndeces function,

Dive deep into the indexes of FFT outputs

Let’s create an artificial series with a hidden cycle element which repeat every 3 days

d_n <- 30
cycle_duration_in_period <- 3
d <- rnorm(d_n)
d[which(0:(d_n-1)%%cycle_duration_in_period==0)] <- 3
plot.ts(d, main="The Seires")

Calculate the FFT of this series and plot

fftd <- fft(d)

x_labels <- paste0(0:(length(d)-1), "/", length(d))
plot(0:(d_n-1), Re(fftd), type = "o", xlab = "Index", ylab = "Re(fftd)", xaxt = "n", main="The FFT Series")
axis(1, at = 0:(d_n-1), labels = x_labels)

We can see the FFT successfully detects the energy of the frequency at \(f = 10/30\), or \(2\pi/3\) radian

Insight from the chart: if we want to create a frequency band covering the frequency component of this 3-day cycle, we want the FFT elements of index 11 ~ 21, corresponding to 10/30 ~ 20/30 (Note: index 1 is 0/30)

The required indexes are:
 [1] 11 12 13 14 15 16 17 18 19 20 21

However, if we apply the method used in (Chatziantoniou, I., Gabauer, D., & Gupta, R., 2021), which implies a specification of below to generate the required FFT indexes

getIndeces(d_n, up = pi, down = pi / cycle_duration_in_period)
 [1]  6  7  8  9 10 11 12 13 14 15 17 18 19 20 21 22 23 24 25 26

which is different from the expected indexes

 [1] 11 12 13 14 15 16 17 18 19 20 21

The difference comes from that the getIndeces function samples the FFT index in the first half of the series, and then further include the mirrored negative frequencies.

The first half of the FFT index contains frequency components for frequency from \(0/N = 0\) to \(\frac{(N/2)}{N}\) = \(\frac{1}{2}\)

Sampling within the range from low frequency \(\pi/K\) to high(est) frequency \(\pi\), correspond to the frequency components that contains FULL CYCLES of frequency from \(\frac{1}{2K}\), (i.e. \(1/6\) if \(K = 3\)) to \(1/2\)

The only sensible way for the Chatziantoniou, et al.’s interpretation to be valid is that they are talking above HALF CYCLES, i.e. the duration from the “PEAK” to the “TROUGH” of the cycle. (without having to going back from TROUGH to PEAK). Only in such interpretation the idea of “1” day / “1-5” days “frequency” (actually “[half-]cycling period”) makes sense in the context of Fourier Transform.

Revisit the concept “Nyquist Frequency”: for the sampling rate of daily data, the highest measurable frequency is 1/2 day, i.e. a 2-day cycle. Any measurement of higher frequency results in “aliasing”, corresponding to the negative frequency of the data.

Notes

Relationship with Standard DFT

\[ \begin{equation}X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2\pi \frac{k}{N} n}\end{equation} \]

  • Rearranging \[ \begin{equation}X_k = \sum_{n=0}^{N-1} e^{-i 2\pi \frac{k}{N} n} x_n \end{equation} \]

  • Match Notation

    • \(n => h\)
    • \(X_n => \Psi_h\)
    • \(2\pi\frac{k}{N} => w\) \[ \begin{equation}X_k = \sum_{h=0}^{N-1} e^{-i w h} \Psi_h \end{equation} \]

Spectral Density

Represent \(x_t\) as \(MA(\infty)\) Filtered Series: \(\begin{equation} \mathbf{x}_t = \sum_{h=0}^{\infty} \Psi_h \boldsymbol{\epsilon}_{t-h} \end{equation}\)

\[ \begin{align*} \mathbf{x}_t &= \sum_{j=0}^{\infty} \Psi_j \boldsymbol{\epsilon}_{t-j} \\ \mathbf{x}_{t-h} &= \sum_{k=0}^{\infty} \Psi_k \boldsymbol{\epsilon}_{t-h-k} \\ \\ \text{So,} \\ \\ \mathbf{x}_{t-h}' &= \sum_{k=0}^{\infty} \boldsymbol{\epsilon}_{t-h-k}' \Psi_k' \end{align*} \]

Therefore, \[ \begin{align*} E\left[\mathbf{x}_t \mathbf{x}_{t-h}'\right] &= E\left[ \left( \sum_{j=0}^{\infty} \Psi_j \boldsymbol{\epsilon}_{t-j} \right) \left( \sum_{k=0}^{\infty} \boldsymbol{\epsilon}_{t-h-k}' \Psi_k' \right) \right] \\ \end{align*} \]

Expand the sums: \[ \begin{align*} &= \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \Psi_j\, E\left[ \boldsymbol{\epsilon}_{t-j} \boldsymbol{\epsilon}_{t-h-k}' \right] \Psi_k' \end{align*} \]

Use the White Noise Properties, the term is non-zero only when \(j = k + h\), \[ \begin{equation} E\left[\boldsymbol{\epsilon}_{t-j} \boldsymbol{\epsilon}_{t-h-k}'\right] = \begin{cases} \Sigma & \text{if } t-j = t-h-k \\ 0 & \text{otherwise} \end{cases} \end{equation} \]

Set \(k = j - h\), \[ \begin{equation} E\left[\mathbf{x}_t \mathbf{x}_{t-h}'\right] = \sum_{j=0}^{\infty} \Psi_j \Sigma \Psi_{j-h}' \end{equation} \]

Substituting to the Spectral Density Formula,

\[ \begin{align*} \sum_{h=-\infty}^{\infty} E(\mathbf{x}_t \mathbf{x}_{t-h}') e^{-i \omega h} = \sum_{h=-\infty}^{\infty} \sum_{j=0}^{\infty} \Psi_j \Sigma \Psi_{j-h}' e^{-i \omega h} \\ = \sum_{j=0}^{\infty} \Psi_j \Sigma \left( \sum_{h=-\infty}^{\infty} \Psi_{j-h}' e^{-i \omega h} \right) \end{align*} \]

Key References

  • BK (2018) Measuring the Frequency Dynamics of Financial Connectedness and Systemic Risk. Journal of Financial Econometrics

Footnotes

  1. Discrete Fourier Transform:↩︎