How to find a good fit for semi-sinusoidal model in R?

Question

I want to assume that the sea surface temperature of the Baltic Sea is the same year after year, and then describe that with a function / linear model. The idea I had was to just input year as a decimal number (or num_months/12) and get out what the temperature should be about that time. Throwing it into lm() function in R, it doesn't recognize sinusoidal data so it just produces a straight line. So I put the sin() function within a I() bracket and tried a few values to manually fit the function, and that gets close to what I want. But the sea is warming up faster in the summer and then cooling off slower in the fall... So the model is wrong the first year, then gets more correct after a couple of years, and then in the future I guess it becomes more and more wrong again.

How can I get R to estimate the model for me, so I don't have to guess numbers myself? The key here is that I want it to produce the same values year after year, not just be correct for one year. If I knew more about math, maybe I could guesstimate it as something like a Poisson or Gaussian instead of sin(), but I don't know how to do that either. Any help to get closer to a good answer would be greatly appreciated.

Here is the data I use, and the code to show results so far:

# SST from Bradtke et al 2010
ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)
SSTlm <- lm(SST$Degrees ~ I(sin(pi*2.07*SST$ToY)))
summary(SSTlm)
plot(SST,xlim=c(0,4),ylim=c(0,17))
par(new=T)
plot(data.frame(ToY=SST$ToY,Degrees=8.4418-6.9431*sin(2.07*pi*SST$ToY)),type="l",xlim=c(0,4),ylim=c(0,17))

Answer 1

It can be done with linear regression -

You just need both a $\sin$ and a $\cos$ term at each frequency.

The reason why you can use a $\sin$ and $\cos$ term in a linear regression to handle seasonality with any amplitude and phase is because of the following trigonometric identity:

A 'general' sine wave with amplitude $A$ and phase $\varphi$, $A \sin (x + \varphi)$, can be written as the linear combination $a\sin x+b\cos x$ where $a$ and $b$ are such that $A=\sqrt{a^2+b^2}$ and $\sin\varphi = \frac{b}{\sqrt{a^2+b^2}}$. Let's see that the two are equivalent:

\begin{eqnarray} a \sin(x) + b \cos(x) &=& \sqrt{a^2+b^2} \left(\frac{a}{\sqrt{a^2+b^2}} \sin(x) + \frac{b}{\sqrt{a^2+b^2}} \cos(x)\right)\\ &=& A\left[\sin(x)\cos(\varphi) + \cos(x)\sin(\varphi)\right]\\ &=& A\sin(x+\varphi)\,\text{.} \end{eqnarray}

Here's the 'basic' model:

 SSTlm <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY),data=SST)
 summary(SSTlm)

[snip]

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)              8.292      0.135   61.41   <2e-16 *** 
sin(2 * pi * ToY)       -5.916      0.191  -30.98   <2e-16 ***  
cos(2 * pi * ToY)       -4.046      0.191  -21.19   <2e-16 *** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9355 on 45 degrees of freedom
Multiple R-squared: 0.969,      Adjusted R-squared: 0.9677 
F-statistic: 704.3 on 2 and 45 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylim=c(1.5,16.5),data=SST)
 lines(SST$ToY,SSTlm$fitted,col=2)

Edit: Important note - the $2\pi\,t$ term works because the period of the function has been set up so that one period = 1 unit of $t$. If the period is different from 1, say the period is $\omega$, then you need $(2\pi/\omega)\, t$ instead.

Here's the model with the second harmonic:

 SSTlm2 <- lm(Degrees ~ sin(2*pi*ToY)+cos(2*pi*ToY)
                        +sin(4*pi*ToY)+cos(4*pi*ToY),data=SST)
 summary(SSTlm2)

[snip]

Coefficients:
                  Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        8.29167    0.02637  314.450  < 2e-16 ***  
sin(2 * pi * ToY) -5.91562    0.03729 -158.634  < 2e-16 ***  
cos(2 * pi * ToY) -4.04632    0.03729 -108.506  < 2e-16 ***  
sin(4 * pi * ToY)  1.21244    0.03729   32.513  < 2e-16 ***  
cos(4 * pi * ToY)  0.33333    0.03729    8.939 2.32e-11 ***  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1827 on 43 degrees of freedom
Multiple R-squared: 0.9989,     Adjusted R-squared: 0.9988 
F-statistic:  9519 on 4 and 43 DF,  p-value: < 2.2e-16 

 plot(Degrees~ToY,ylab="Degrees",xlab="ToY",ylim=c(1.5,16.5),data=SST)
 lines(SSTlm2$fitted~ToY,col=2,data=SST)

... and so forth, with 6*pi*ToY etc. If there was a tiny bit of noise in the data I'd probably stop with this second model though.

With enough terms, you can exactly fit asymmetric and even jagged periodic sequences, but the resulting fits may 'wiggle'. Here's an asymmetric function (it's a sawtooth -) added to a scaled version of your periodic function), with third (red) and fourth (green) harmonics. The green fit is on average a little closer but "wiggly" (even when the fit goes through every point, the fit may be very wiggly between points).

The periodicity here means there's only 12 d.f. available for a seasonal model in the data. With the intercept in the model, you only have enough degrees of freedom for 11 additional seasonal parameters. Since you are adding two terms with each harmonic, the last harmonic you can fit will only allow you one of them for the last term, the sixth harmonic (and that one has to be a $\cos$; the $\sin$ term will be all-zero, while the cos alternates between 1 and -1).

If you want fits that are smoother than this approach produces on non-smooth series, you may want to look into periodic spline fits.

Yet another approach is to use seasonal dummies, but the sin/cos approach is often better if it's a smooth periodic function.

This kind of approach to seasonality can also adapt to situations where seasonality is changing, such as using trigonometric or dummy seasonality with state-space models.

---

While the linear model approach discussed here is simple to use, one advantage of @COOLSerdash's nonlinear regression approach is that it can deal with a much wider range of situations - you don't have to change much before you're in a situation where linear regression is no longer suitable but nonlinear least-squares may still be used (having an unknown period would be one such case).

Answer 2

The temperature you provide in your question repeats exactly every year. I suspect this aren't really measured temperatures over four years. In your example, you wouldn't need a model, because the temperatures just repeat exactly. But otherwise you could use the nls function to fit a sine curve:

ToY <- c(1/12,2/12,3/12,4/12,5/12,6/12,7/12,8/12,9/12,10/12,11/12,12/12,13/12,14/12,15/12,16/12,17/12,18/12,19/12,20/12,21/12,22/12,23/12,24/12,25/12,26/12,27/12,28/12,29/12,30/12,31/12,32/12,33/12,34/12,35/12,36/12,37/12,38/12,39/12,40/12,41/12,42/12,43/12,44/12,45/12,46/12,47/12,48/12)
Degrees <- c(3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5,3,2,2.2,4,7.6,13,16,16.1,14,10.1,7,4.5)
SST <- data.frame(ToY, Degrees)

par(cex=1.5, bg="white")
plot(Degrees~ToY,xlim=c(0,4),ylim=c(0,17), pch=16, las=1)

nls.mod <-nls(Degrees ~ a + b*sin(2*pi*c*ToY), start=list(a = 1, b = 1, c=1))

co <- coef(nls.mod) 
f <- function(x, a, b, c) {a + b*sin(2*pi*c*x) }

curve(f(x, a=co["a"], b=co["b"], c=co["c"]), add=TRUE ,lwd=2, col="steelblue")

But the fit isn't very good, especially at the beginning. It seems that your data cannot be adequately modelled by a simple sine curve. Maybe a more complex trigonometrical function will do the trick?

nls.mod2 <-nls(Degrees ~ a + b*sin(2*pi*c*ToY)+d*cos(2*pi*e*ToY), start=list(a = 1, b = 1, c=1, d=1, e=1))

co2 <- coef(nls.mod2) 
f <- function(x, a, b, c, d, e) {a + b*sin(2*pi*c*x)+d*cos(2*pi*e*x) }

curve(f(x, a=co2["a"], b=co2["b"], c=co2["c"], d=co2["d"], e=co2["e"]), add=TRUE ,lwd=2, col="red")

The red curve fits the data better. With the nls function, you can put in the model that you think is appropriate.

Or maybe you could make use the the forecast package. In the example below, I have assumed that the time series started in January 2010:

library(forecast)

Degrees.ts <- ts(Degrees, start=c(2010,1), frequency=12)

Degree.trend <- auto.arima(Degrees.ts)

degrees.forecast <- forecast(Degree.trend, h=12, level=c(80,95), fan=F)

plot(degrees.forecast, las=1, main="", xlab="Time", ylab="Degrees")

Because the data is deterministic, no confidence bands are shown.