The following are the issues that I summarize succinctly in sequence, as each of them are an interesting problem in itself and then follow it up with my solutions:
i) How to unlag two time-series so as to lose the least data due to the lagging, and at the same time to maximize the cross-correlation or spectral coherence.
ii) how would this be done when you are dealing with multivariate time-series. As it's greater than the two, its like solving a rubiks cube optimally while lagging various time-series in the dataset so as to preserve a summed up measure of (i).
iii) How to compute a weighted-coherence?
I will skip ii) and answer i) and iii) due to reasons of conflict of interest and research collaboration/confidential agreements on my work on ii). I understand that iii) happens to be the central question to start with followed by i).
My solution for iii)
Pre-Processing Steps: Say you had two time -series t1 and t2. Determine a CCF maximizing lag, say optLag and unlag t1 and t2 as unlaggedTS = lag(t2,optLag) and perform a time-series union of the dataset as unionTS = ts.union(t1,unlaggedTS).
Following that, as the time-series unlagging would produce a few NA's as you lose entries by lagging, you need to set these values to zero in unionTS as required-meaning that nothing was observed at the NA points.
Step 2:
Now, obtain an abinded time-series, with the span parameters as required.
bindTSPair= abind(((spectrum(unionTS,span=c(16,16),plot=F)$spec)[,1]),((spectrum(unionTS,span=c(16,16),plot=F)$spec)[,2]))
Following it, just normalize the data in bindTSPair as say:
normalizedTS=round(data.Normalization(c(bindTSPair),type="n4"),6)
Then obtained the following two halves:
TSHalf1=normalizedTS[1:(length(bindTSPair)/2)]
TSHalf2=normalizedTS[((length(bindTSPair)/2)+1):(length(bindTSPair))]
and then finally obtain 'a weighted' version of spectral coherence as follows:
weightedCoherence=sum((10*(TSHalf1 + TSHalf2)) * (spectrum(unionTS,span=c(10,10),plot=F)$coh))
You may change the weighting as required or using 'prior' info.
My solution for i) to obtain a CCF maximizing lag:
for(i in 2:ncol(myTS))
{
for(j in 2:ncol(myTS))
{
r=ccf(myTS[,i],myTS[,j],lag.max=1000,plot=TRUE)
tp=sort((r$acf),index.return=TRUE)$ix
tn=sort((r$acf),index.return=TRUE,decreasing=T)$ix
cn=r$acf[tn[length(tn)]]
cp=r$acf[tp[length(tp)]]
if(abs(cp) > abs(cn))
{
temp=r$lag[tp[length(tp)]]
tcc=cp
}
if(abs(cn) > abs(cp))
{
temp=r$lag[tn[length(tn)]]
tp=tn
tcc=cn
}
if(abs(cn) > abs(cp))
{
temp=r$lag[tp[length(tp)]]
}
show(i)
if(length(tp) != 0)
{
mlag[i,j]=temp
mcc[i,j]=tcc
}
}
}
mlag=mlag[2:ncol(myTS),2:ncol(myTS)]
mlag=data.frame(mlag)
It all boils down on how would you want to process a time series as to breakdown its components as to use these for later prediction or classification.
For one ARIMA is a parametric method (assumption of a fixed distribution) modeling a stationary time series based on static ARMA terms while with wavelets you model a wavelet function by selecting a list of characterictis you want the function to have as to best approximate a signal (wavelets can model a non stationary as well as stationary). In wavelets the length of the filter the number of vanishing moments and the symetry of the mother wavelet vs the signal will define how good the function is in providing representation of time resolution and frequency resolution. (existing a trade off between this two terms as you get shorter filters you get better in time resolution while decrease the frequency resolution.
In ARIMA you aproximate the signal components by selecting ARMA terms from acf and pacf while for wavelets you aproximate the signal by selecting a mother wavelet, filter length and number of vanishing moments. (the number of decomposition levels provides capturing certain frequencies per level). So an stationary ARMA model would have its pseudo equivalent for a non stationary signal by performing a wavelet transform only at one scale obtaining scaling and detail coefficients (sort of equivalent to having an ARMA with dynamic terms). The results for the two methods will be different, for as in Arima you have static ARMA terms across time, in wavelet decomposition you will have frequency ranges across time which will not correspond to the ARMA terms, although you could calculate what are the frequencies that the ARMA terms represent and approximate these by using wavelets, the dynamic approach of the wavelet will ensure that there is no one specific frequency across time where the ARMA terms will be represented).
In modeling a time series the first issue is to define the characteristics of the signal and the second is in knowing in what part of the signal we are most interested in. (modeling the entire signal or extracting a portion of it).
Which modeling method will provide a better representation of the signal or portions of it you are interested in, depends on the ability of the modeler, the properties of the signal (how easy is to model it) Sometimes an stockastic model such us ARMA or other modeling method model will be sufficient, other times using a built in transform by selecting carefully its terms will give better representation.
Best Answer
The Fourier transform of discrete data with sampling rate $\Delta t = 1$ is periodic with period $2\pi$ (if we're talking about angular frequency, $\omega$) OR with period $1$ (if we're talking about linear frequency, $f$). The relationship between angular and linear frequency is $\omega = 2\pi f$. We'll assume for simplicity that $\Delta t = 1$ for now and I'll be using linear frequency, $f$.
In fact, with real-valued data the Fourier transform (which is complex) is conjugate-symmetric about $\omega = f = 0$. If $X(f)$ is the Fourier transform of $x(t) \in \mathbb{R}$, then $X(f) = X^{*}(-f)$. The spectrum is proportional to $|X(f)|^{2}$, so the spectrum is symmetric about $f = 0$.
So, with a discretely sampled data ("index limited") we end up with a "band limited" Fourier transform. The Nyquist frequency (the highest frequency which is not an alias of lower frequencies. The below is an example of aliasing (can't tell one sinusoid from another - see Aliasing on Wikipedia).
If $\Delta t \neq 1$ then the Nyquist frequency is $1 / (2\Delta t)$ and the frequencies where the spectrum is calculated (if you have $N$ data points) is
Usually you want to zeropad for the FFT, so you add some zeroes on at the end.
One of the better options for estimating the spectrum is to use the Multitaper Method (package multitaper in R from CRAN).
Sorry, this was all over the place... An example is below. The plot that it creates demonstrates the bias (broadband bias) properties of the periodogram vs. the multitaper.
An example I thoroughly enjoy comes from Percival and Walden, page 46 - an order 4 autoregressive simulation (so we know the theoretical spectrum): $$x_{t} = 2.7607x_{t-1} - 3.8106x_{t-2} + 2.6535x_{t-3} - 0.9238x_{t-4}$$