I am trying to remove seasonality from data. I tried the non-linear trend using the code: trend=lm(NH3cH6~t+cos.t+sin.t). The plot was shown as following:

However, as you can see, the second peak of the trend does not quite fit the data. I want to shift the second peak of the trend a little bit, but I don't know how to modify the code to achieve this. Can anyone help me to improve the fitting of the trend?

• You've fitted a model with period $2\pi$. What's the period of your data? And what's the frequency of the observations? Commented Nov 20, 2013 at 11:27
• Perhaps too obvious, but other small problems too. Fitted function if linear plus sinusoid should be smoother than you have shown and be defined for the whole data range. Commented Nov 20, 2013 at 13:52
• They are two-year daily data. Each peak indicates one year. Thanks! @RobHyndman Commented Nov 20, 2013 at 19:20
• Chemical comment: the "C" in NH$_3$CH$_6$ presumably means carbon and should be upper case. Commented Nov 21, 2013 at 9:42

Your model fits a linear time trend plus a first order Fourier series approximation for the seasonality. Since you don't define cos.t or sin.t, it is not possible to tell what seasonal period you have used. You can fit a model with correct seasonal period using

t <- 1:length(NH3cH6)
cos.t <- cos(2*pi*t/365)
sin.t <- sin(2*pi*t/365)
trend <- lm(NH3cH6 ~ t + cos.t + sin.t)


However, you may need a higher order Fourier series approximation, and there are facilities in the forecast package in R for doing that.

First, make sure the data is a time series object of the correct frequency:

NH3cH6 <- ts(NH3cH6, frequency=365)


Generate the Fourier terms for the regression (using 3 terms here, but choose the best number by minimizing the AIC of the fitted regression model).

library(forecast)
X <- fourier(NH3cH6,3)


Fit the regression model

fit  <- tslm(NH3cH6 ~ X + trend)


The tslm command will understand what trend means, and will make sure the residuals and fitted values are ts objects with time series characteristics.

To see the result:

plot(NH3cH6)
lines(fitted(fit), col="red")


I'm not at all sure that a trigonometric wave will be a good fit for these data. Look at how the peaks seem so narrow and pointy and the troughs seem so wide and flat. It's also not clear to me why you would include the sum of the sin and the cos functions. Do you have some theoretical reason behind the formula that you started with?

If you wish to try fitting a simple trig function, you might want to explore nonlinear regression instead. For example,

# a is the vertical shift (height of horizontal center line)
# b is the amplitude (half the vertical distance from trough to peak)
# c is the phase shift (distance from first trough to x=0)
# d is the period (distance from peak to peak)
trend <- nls(NH3cH6 ~ a + b * sin(2*pi/d*(c+t)), start=list(a=60, b=60, c=-80, d=400))


Alternatively, you could fit a smooth line to the data, e.g., using the R functions loess.smooth() or gam() from the mgcv package.

• A trigonometric model can fit periodic series that have narrow peaks and wide troughs ... if you use more than one harmonic. You use sums of sin and cos functions because the amplitude and phase are unknown; the coefficients on sin and cos are related to the amplitude and phase. The benefit of working with sin and cos is that you end up with linear rather than nonlinear models. Both these issues are discussed here Commented Nov 21, 2013 at 7:40
• In addition, using a root or log link and a generalized linear model might be a reasonable way to cope with the mild asymmetry you identify, which seems more a consequence of skewness of distribution than anything else. Commented Nov 21, 2013 at 9:41
• Thanks @Glen_b for that link. Very interesting. And a useful approach, as long as the period is known. Commented Nov 21, 2013 at 14:11
• @JeanV.Adams As mentioned at the end of that answer, there's discussion of unknown period here Commented Nov 21, 2013 at 23:06

Fitting trigonometric models seems nice because you can describe relatively complicated trends with relatively few parameters. However, the lessons learned are always the same: don't do it! Even as a pure mathematics problem, there are plenty of convergence issues involved, identifying correct periodicity, amplitude, etc. Add the randomness of a statistical experiment to the equation and you have a mess on your hands.

Smoothers are great at accounting for unknown functions of time, trigonometric or otherwise. The basic boxcar smoother is intuitive and easy to apply in this case, but you can also consider smoothing splines or a LOESS smooth. With a sufficient amount of data, it's possible to fit effects at an annual, quarterly, or monthly level which are consistently estimated over a long period of time.

I am vaguely aware that in advanced atmospheric sciences, the Kalman Filter is somehow applied to estimate more granular seasonal trends. Caveat emptor! But here is an applied paper which looks at just this problem. My basic understanding indicates that the Kalman filter is well suited to this kind of problem since it jointly estimates state change and steady state effects in a model. Please no engineers critique my explanation. The benefit of the filter over smoothers is that they are used to predict values out in the future (or "forecasting" as it's called in such fields). As a result, the Kalman estimated trends are believably more "consistent" whereas smoothers are not really meant to infer many differences for specific values (though you can do it that way, it's not a conventional approach and does have shortcomings, aka "Curse of Dimensionality").