Kolmogorov-Zurbenko third-order periodogram and smoothing method

1 2 3 | ```
kztp(x,m,k,p=1,n=1,rp1=0,rp2=0.5,cp1=0,cp2=0.5)
variation.kztp(pg, K=dim(pg)[1])
smooth.kztp(pg,c,K=dim(pg)[1])
``` |

`x` |
The raw time series |

`m` |
The length of the window size for a regular Fourier transform |

`k` |
The number of iterations for the KZFT |

`p` |
The distance between two successive intervals as a percentage of the total length of the time series |

`n` |
The sampling frequency rate as a multiplication of the Fourier frequencies |

`rp1` |
The left bound of the frequencies at which the third-order periodogram will be calculated in one frequency dimension |

`rp2` |
The right bound of the frequencies at which the third-order periodogram will be calculated in one frequency dimension |

`cp1` |
The left bound of the frequencies at which the third-order periodogram will be calculated in the other frequency dimension |

`cp2` |
The right bound of the frequencies at which the third-order periodogram will be calculated in the other frequency dimension |

`pg` |
The raw third-order periodogram of a time series |

`K` |
Half of the maximum bandwidth of the spectral window being allowed |

`c` |
A prespecified percentage of total variation |

Kolmogorov-Zurbenko third-order periodogram is calculated based on the Kolmogorov-Zurbenko Fourier transform. The smoothing algorithm is very similar to the DZ algorithm used in smoothing the periodogram. It is also an adaptive algorithm which allows the bandwidth of the spectral window to vary according to the smoothness of the underlying bispectrum density. The spectral window being used here is a square window. In the algorithm, the bandwith of the square window is extended until the squared variation of the third-order periodogram within the window reaches a prespecified percentage of the total variation of the raw third-order periodogram.

I. G. Zurbenko, 1986: The spectral Analysis of Time Series. North-Holland, 248 pp.

I. G. Zurbenko, P. S. Porter, Construction of high-resolution wavelets, Signal Processing 65: 315-327, 1998.

W. Yang, I. G. Zurbenko, A semi-adaptive smoothing algorithm in bispectrum estimation, Proceedings of the American Statistical Association, Seattle, 2006.

1 2 3 4 5 6 7 8 9 10 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 56 57 58 59 60 | ```
#example 1
t<-1:10000
y<-sin(2*pi*0.1*t)+sin(2*pi*0.2*t)+rnorm(10000,0,3)
pg3.y<-kztp(y,100,1,0.5,1)$mod
pg3.y<-log(pg3.y)
spg3.y<-smooth.kztp(pg3.y,0.01,20)$bispectrum
omega<-seq(0,1,length=101)[2:51]
par(mfrow=c(2,1))
persp(omega, omega, pg3.y, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed", xlab = "Frequency (cycles/time unit)",
ylab = "Frequency (cycles/time unit)", main = "3rd-order Periodogram")
persp(omega, omega, spg3.y, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed", xlab = "Frequency (cycles/time unit)",
ylab = "Frequency (cycles/time unit)", main = "Smoothed 3rd-order Periodogram")
filled.contour(omega,omega,pg3.y,xlab="Frequency (cycles/time unit)",ylab="Frequency (cycles/time unit)",
main="3rd-order Periodogram")
filled.contour(omega,omega,spg3.y,xlab="Frequency (cycles/time unit)",ylab="Frequency (cycles/time unit)",
main="Smoothed 3rd-order Periodogram3")
#example 2
#effect of KZFT
t<-1:50000
y<-1.1*sin(2*pi*0.0339*t)+7*sin(2*pi*0.0366*t)+2*rnorm(50000,0,1)
pg3.y1<-kztp(y,1000,1,0.5,1,rp1=0.02,rp2=0.05,cp1=0.02,cp2=0.05)$mod
pg3.y2<-kztp(y,1000,3,0.5,1,rp1=0.02,rp2=0.05,cp1=0.02,cp2=0.05)$mod
pg3.y1<-log(pg3.y1)
pg3.y2<-log(pg3.y2)
spg3.y1<-smooth.kztp(pg3.y1,0.01,10)$bispectrum
spg3.y2<-smooth.kztp(pg3.y2,0.01,10)$bispectrum
omega<-seq(0,1,length=1001)[21:51]
par(mfrow=c(2,2))
persp(omega, omega, pg3.y1, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed", xlab = "Frequency (cycles/time unit)",
ylab = "Frequency (cycles/time unit)", main = "Original, m=1000, k=1")
persp(omega, omega, pg3.y2, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed", xlab = "Frequency (cycles/time unit)",
ylab = "Frequency (cycles/time unit)", main = "Original, m=1000, k=3")
persp(omega, omega, spg3.y1, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed", xlab = "Frequency (cycles/time unit)",
ylab = "Frequency (cycles/time unit)", main = "Smoothed, m=1000, k=1")
persp(omega, omega, spg3.y2, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
ltheta = 120, shade = 0.75, ticktype = "detailed", xlab = "Frequency (cycles/time unit)",
ylab = "Frequency (cycles/time unit)", main = "Smoothed, m=1000, k=3")
filled.contour(omega,omega,pg3.y1,xlab="Frequency (cycles/time unit)",ylab="Frequency (cycles/time unit)",
main="Original, m=1000, k=1")
filled.contour(omega,omega,pg3.y2,xlab="Frequency (cycles/time unit)",ylab="Frequency (cycles/time unit)",
main="Original, m=1000, k=3")
filled.contour(omega,omega,spg3.y1,xlab="Frequency (cycles/time unit)",ylab="Frequency (cycles/time unit)",
main="Smoothed, m=1000, k=1")
filled.contour(omega,omega,spg3.y2,xlab="Frequency (cycles/time unit)",ylab="Frequency (cycles/time unit)",
main="Smoothed, m=1000, k=3")
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.