#' 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))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.