R/DFA_code.R

In wiaidp/MDFA: Signalextraction and Forecasting Based on the Multivariate Direct Filter Approach (MDFA)

#' Simple function for computing DFT and periodogram (based on univariate DFA approach): cannot handle multivariate signal extraction
#'
#' @param x Data (time series)
#' @param plot_T Boolean (generate a plot of the periodogram: T/F)
#' @return DFT Discrete Fourier Transform
#' @return per Periodogram (squared absolute DFT)
#' @export
#' @import graphics



per<-function(x,plot_T)
{
  len<-length(x)
  per<-0:(len/2)
  DFT<-per

  for (k in 0:(len/2))
  {
    cexp <- exp(1.i*(1:len)*2*pi*k/len)
    DFT[k+1]<-sum(cexp*x*sqrt(1/(2*pi*len)))
  }
# Frequency zero receives weight 1/sqrt(2)
#   The periodogram in frequency zero appears once only whereas all other frequencies are doubled

# This is omitted now in order to comply with MDFA
#   We now change the periodogram in the dfa estimation routines
#  DFT[1]<-DFT[1]/sqrt(2)
# Weighths wk: if length of data sample is even then DFT in frequency pi is scaled by 1/sqrt(2) (Periodogram in pi is weighted by 1/2)
  if (abs(as.integer(len/2)-len/2)<0.1)
    DFT[k+1]<-DFT[k+1]/sqrt(2)
  per<-abs(DFT)^2
  if (plot_T)
  {
    par(mfrow=c(2,1))
    plot(per,type="l",axes=F,xlab="Frequency",ylab="Periodogram",
         main="Periodogram")
    axis(1,at=1+0:6*len/12,labels=c("0","pi/6","2pi/6","3pi/6",
                                    "4pi/6","5pi/6","pi"))
    axis(2)
    box()
    plot(log(per),type="l",axes=F,xlab="Frequency",ylab="Log-periodogram",
         main="Log-periodogram")
    axis(1,at=1+0:6*len/12,labels=c("0","pi/6","2pi/6","3pi/6",
                                    "4pi/6","5pi/6","pi"))
    axis(2)
    box()
  }
  return(list(DFT=DFT,per=per))
}








#' This function computes MSE (mean-square) DFA-solutions: no customization, univariate problems only, no regularization, no constraints
#'
#' @param L Filter length
#' @param periodogram Periodogram (as generated by function per)
#' @param Lag Backcast (Lag>0), Nowcast (Lag=0) or Forecast (Lag<0)
#' @param Gamma Generic target specification: typically symmetric Lowpass (trend) or Bandpass (cycle) filters. Highpass and anticipative allpass (forecast) can be specified too.
#' @return b optimal (MSE) filter coefficients
#' @return trffkt Complex transfer function of optimal (univariate MSE) filter
#' @export
#'
dfa_ms<-function(L,periodogram,Lag,Gamma)
{
# Frequency 0 appears once only whereas all other frequencies are doubles: therefore we halve periodogram in frequency zero
  periodogram[1]<-periodogram[1]/2
  K<-length(periodogram)-1
  X<-exp(-1.i*Lag*pi*(0:(K))/(K))*rep(1,K+1)*sqrt(periodogram)
  X_y<-exp(-1.i*Lag*pi*(0:(K))/(K))*rep(1,K+1)
  for (l in 2:L)          #l<-L<-21
  {
    X<-cbind(X,(cos((l-1-Lag)*pi*(0:(K))/(K))+
                  1.i*sin((l-1-Lag)*pi*(0:(K))/(K)))*sqrt(periodogram))
    X_y<-cbind(X_y,(cos((l-1-Lag)*pi*(0:(K))/(K))+
                      1.i*sin((l-1-Lag)*pi*(0:(K))/(K))))
  }
  xtx<-t(Re(X))%*%Re(X)+t(Im(X))%*%Im(X)
  # MA-Filtercoefficients
  b<-as.vector(solve(xtx)%*%(t(Re(X_y))%*%(Gamma*periodogram)))
  # Transferfunction
  trffkt<-1:(K+1)
  trffkt[1]<-sum(b)
  for (k in 1:(K))#k<-1
  {
    trffkt[k+1]<-(b%*%exp(1.i*k*(0:(length(b)-1))*pi/(K)))
  }
  return(list(b=b,trffkt=trffkt))
}







#' This function computes MSE (mean-square) DFA-solutions as well as customized filters: it replicates the function dfa_ms if lambda=eta=0 and if i1=i2=F. It can solve univariate signal extraction problems only. It can handle constraints. Regularization is not possible.
#'
#' @param L Filter length
#' @param lambda Customization parameter: Timeliness is emphasized in the ATS-trilemma if lambda>0
#' @param periodogram Periodogram (as generated by function per)
#' @param Lag Backcast (Lag>0), Nowcast (Lag=0) or Forecast (Lag<0)
#' @param Gamma Generic target specification: typically symmetric Lowpass (trend) or Bandpass (cycle) filters. Highpass and anticipative allpass (forecast) can be specified too.
#' @param eta Customization parameter: Smoothness is emphasized in the ATS-trilemma if eta>0
#' @param cutoff Specifies start-frequency in stopband from which Smoothness is emphasized (corresponds typically to the cutoff of the lowpass target). Is not used if eta=0.
#' @param i1 Boolean. If T a first-order filter constraint in frequency zero is obtained: amplitude of real-time filter must match weight_constraint (handles integration order one)
#' @param i2 Boolean. If T a second-order filter constraint in frequency zero is obtained: time-shift of real-time filter must match target (together with i1 handles integration order two)
#' @return b optimal filter coefficients: MSE-design if lambda=eta=0
#' @return trffkt Complex transfer function of optimal (univariate) filter
#' @export
#'
dfa_analytic<-function(L,lambda,periodogram,Lag,Gamma,eta,cutoff,i1,i2)
{
# Frequency 0 appears once only whereas all other frequencies are doubles: therefore we halve periodogram in frequency zero
  periodogram[1]<-periodogram[1]/2
  # Impose meaningful parameter restrictions
  lambda<-abs(lambda)
  eta<-abs(eta)
  K<-length(periodogram)-1
  # Define the amplitude weighting function weight_h (W(omega_k))
  omega_Gamma<-as.integer(cutoff*K/pi)
  if ((K-omega_Gamma+1)>0)
  {
    weight_h<-periodogram*(c(rep(1,omega_Gamma),(1:(K-omega_Gamma+1))^(eta)))
  } else
  {
    weight_h<-periodogram*rep(1,K+1)
  }
  # First order filter restriction: assigning a `large' weight to frequency zero
  if (i1)
    weight_h[1]<-max(1.e+10,weight_h[1])

  X<-exp(-1.i*Lag*pi*(0:(K))/(K))*rep(1,K+1)*sqrt(weight_h)
  X_y<-exp(-1.i*Lag*pi*(0:(K))/(K))*rep(1,K+1)
  if (i2)
  {
    # Second order restriction: time shift in frequency zero vanishes
    for (l in 2:(L-1))
    {
      X<-cbind(X,(cos((l-1-Lag)*pi*(0:(K))/(K))-((l-1)/(L-1))*
                    cos((L-1-Lag)*pi*(0:(K))/(K))+
                    sqrt(1+Gamma*lambda)*1.i*(sin((l-1-Lag)*pi*(0:(K))/(K))-((l-1)/(L-1))*
                                                sin((L-1-Lag)*pi*(0:(K))/(K))))*sqrt(weight_h))
      X_y<-cbind(X_y,(cos((l-1-Lag)*pi*(0:(K))/(K))-((l-1)/(L-1))*
                        cos((L-1-Lag)*pi*(0:(K))/(K))+
                        1.i*(sin((l-1-Lag)*pi*(0:(K))/(K))-((l-1)/(L-1))*sin((L-1-Lag)*pi*(0:(K))/(K)))))
    }
    xtx<-t(Re(X))%*%Re(X)+t(Im(X))%*%Im(X)
    # MA-Filterweights
    b<-as.vector(solve(xtx)%*%(t(Re(X_y))%*%(Gamma*weight_h)))
    # the last weight is a function of the previous ones through the second order restriction
    b<-c(b,-sum(b*(0:(length(b)-1)))/(length(b)))
  } else
  {
    for (l in 2:L)
    {
      X<-cbind(X,(cos((l-1-Lag)*pi*(0:(K))/(K))+
                    sqrt(1+Gamma*lambda)*1.i*sin((l-1-Lag)*pi*(0:(K))/(K)))*sqrt(weight_h))
      X_y<-cbind(X_y,(cos((l-1-Lag)*pi*(0:(K))/(K))+
                        1.i*sin((l-1-Lag)*pi*(0:(K))/(K))))
    }
    xtx<-t(Re(X))%*%Re(X)+t(Im(X))%*%Im(X)
    # MA-Filterweights
    b<-as.vector(solve(xtx)%*%(t(Re(X_y))%*%(Gamma*weight_h)))
  }
  # Transferfunction
  trffkt<-1:(K+1)
  trffkt[1]<-sum(b)
  for (k in 1:(K))#k<-1
  {
    trffkt[k+1]<-(b%*%exp(1.i*k*(0:(length(b)-1))*pi/(K)))
  }
  return(list(b=b,trffkt=trffkt))
}