haRmonics: Harmonic analysis for time series

Description

Fits harmonic regression models, that is, computes amplitudes and phase angles in the typical harmonic regression framework. When method=harmR the ordinary least squares method is used, when method=wls_harmR then, weighted least squares are employed. Based on these estimates a harmonic regression function is fitted. Also fits hants, a popular iterative algorithm that computes amplitudes and phase angles in the harmonic regression framework. As part of the iterative algorithm, observations are being excluded from the design matrix of the regression model if the distance between them and the fitted curve exceeds the value of the parameter fitErrorTol. hants is based on implementations with the same name written in Fortran and Matlab computer languages.

Usage

haRmonics(
  y,
  method = c("harmR", "wls_harmR", "hants"),
  sigma = NULL,
  ts = 1:length(y),
  lenPeriod = length(y),
  numFreq,
  HiLo = c("Hi", "Lo"),
  low,
  high,
  fitErrorTol,
  degreeOverDeter,
  delta
)

Value

A list containing:

a.coef: a numeric vector with estimates of cosine coefficients
b.coef: a numeric vector with estimates of sine coefficients
amplitude: a numeric vector with amplitude estimates.
phase: a numeric vector with phase estimates.
fitted: a numeric vector with fitted values via harmonic regression.

Arguments

y: numeric vector containing time series on which harmonic regression will be fitted. Missing values are not allowed.
method: character specifying algorithm to apply: harmR (default), wls_harmR (heteroscedastic model) or hants.
sigma: numeric vector of length lenPeriod containing variance estimates. Default set NULL.
ts: numeric vector of length(y) with the sampling points for y. Default is $t s [i] = i, i = 1, \dots, \code l e n g t h (y)$ .
lenPeriod: numeric giving the length of the base period, reported in samples, e.g. days, dekads, months, years, etc.
numFreq: numeric indicating the total number of frequencies to be used in harmonic regression. For technical reasons, 2*numFreq+1 must be lesser than length(y).
HiLo: character indicating whether high or low outliers must be rejected when method=hants.
low: numeric giving minimum valid value of fitted harmonic regression function when method=hants.
high: numeric giving maximum valid value of fitted harmonic regression function when method=hants.
fitErrorTol: numeric giving maximum allowed distance between observations and fitted curve; if difference between a given observation and its fitted value exceeds fitErrorTol then this observation will not be included in the fitting procedure in the next iteration of the algorithm.
degreeOverDeter: numeric; iteration stops when number of observations equals number of observations for curve fitting plus degreeOverDeter; the latter in turns is by definition length(y) minus min(2 * numFreq+1, length(y)).
delta: numeric (positive) giving a (small) regularization parameter to prevent non-invertible hat matrix (see Details), probably caused by high amplitudes.

Details

Methods harmR and wls_harmR do not allow missing values and utilize parameters y, lenPeriod, numFreq and delta only.

Method hants utilizes all the parameters presented above. This method does not allow missing values. Missing values in y must be substituted by values considerably out of observations range.

References

Roerink, G.J., Menenti, M., Verhoef, W. (2000). Reconstructing cloudfree NDVI composites using Fourier analysis of time series, Int. J. Remote Sensing, 21(9), 1911--1917.

Jakubauskas, M., Legates, D., Kastens, J. (2001). Harmonic analysis of time-series AVHRR NDVI data, Photogrammetric Engineering and Remote Sensing, 67(4), 461--470.

The Matlab implementation of HANTS can be found here.

Examples

Run this code

y <- c(5, 2, 5, 10, 12, 18, 20, 23, 27, 30, 40, 60, 66,
70, 90, 120, 160, 190, 105, 210, 104, 200, 90, 170,
50, 120, 80, 60, 50, 40, 30, 28, 24, 20, 15, 10)
# -------------------------------------------------------------------------- 
fit_harmR <- haRmonics(y = y, numFreq = 3, delta = 0.1)
fitLow_hants <- haRmonics(y = y, method = "hants", numFreq = 3, HiLo = "Lo", 
                         low = 0, high = 255, fitErrorTol = 5, degreeOverDeter = 1, 
                         delta = 0.1)
fitHigh_hants <- haRmonics(y = y, method = "hants", numFreq = 3, HiLo = "Hi", 
                          low = 0, high = 255, fitErrorTol = 5, degreeOverDeter = 1, 
                          delta = 0.1)
plot(y, pch = 16, main = "haRmonics fitting")
lines(fit_harmR$fitted ,lty = 4, col = "green")
lines(fitLow_hants$fitted, lty = 4, col = "red")
lines(fitHigh_hants$fitted, lty = 2, col = "blue")
# -------------------------------------------------------------------------- 
# Substituting missing value by a number outside observations range
# -------------------------------------------------------------------------- 
y1 <- y
y1[20] <- -10

fitLow_hants_missing <- haRmonics(y = y1, method = "hants", numFreq = 3, HiLo = "Lo", 
                                 low = 0, high = 255, fitErrorTol = 5, degreeOverDeter = 1, 
                                 delta = 0.1)
fitHigh_hants_missing <- haRmonics(y = y1, method = "hants", numFreq = 3, HiLo = "Hi", 
                                  low = 0, high = 255, fitErrorTol = 5, degreeOverDeter = 1, 
                                  delta = 0.1)
fit_harmR_missing <- haRmonics(y = y1, numFreq = 3, delta = 0.1)

plot(y1, pch = 16, main = "haRmonics fitting (missing values)", ylim = c(-1,210))
lines(fitLow_hants_missing$fitted, lty = 4, col = "red")
lines(fitHigh_hants_missing$fitted, lty = 2, col = "blue")
lines(fit_harmR_missing$fitted, lty = 4, col = "green")

Run the code above in your browser using DataLab