haRmonics {geoTS}

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
)

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 ts[i] = i, i=1,\ldots, \code{length(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.

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.

Note

lenBasePeriod was used until version 0.1.3, this argument has been replaced by lenPeriod.

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

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")