Why is my sinusodial fit not working well for this data?

Question

I'm trying to find the best fit line for this data below but no matter what I try, the fit line seems to never be able to account for the lower values as shown below.

The x-values are just dates from 1/1/2014 to 7/20/2015 (566 values), but I don't know how to give you guys the y-values. I have it in my Environment but I don't know how to give you that without copying and pasting from the Console output.

This is the code that I'm using to get that fit line:

wb.loglik=function(theta,y,x,null=NA)
{
  a=theta[1]
  b=theta[2]
  c=theta[3]
  d=theta[4]

  if(!is.na(null))
  {
    d=null
  }

  s2=theta[5]
  n=length(y)
  return((-n/2)*log(s2)-1/(2*s2)*sum((y-(a+b*cos(2*pi*((x-c)/d))))^2))
}
result=optim(par=c(mean(wbbcf),sd(wbbcf),1,365.25,var(wbbcf)/2),
fn=wb.loglik,x=Time,y=wbbcf,control=list(fnscale=-1))
theta=result$par
theta
value=result$value
value

This is the code to get the plot above:

plot(date,wbbcf,xlim=c(as.Date("2014-01-01"),as.Date("2015-07-
    20")),ylim=range(c(-3.2,0)),xlab="Date (1/1/2014 to
    7/20/2015)",ylab="Total Net with Storage (bcf)",main="Total Burn with
    Model")
par(new=T)
curve(-0.9740582-0.7857229*cos(2*pi*(x-5.9582996)/385.1581090),1,566,
    ylim=range(c(-3.2,0)),col="blue",xlab="",ylab="",xaxt='n',yaxt='n')

What else can I do to generate a better fit line for this data? Also, if there's a good way to predict future data, I would appreciate help with that as well.

Sorry in advance, if I'm not giving enough information. Please feel free to ask for any information you need and I will promptly edit the post.

EDIT: I have added the work I did with ARIMA below.

I inputted the following code and got the following results:

forecast::auto.arima(wbbcf)
fit.arima = arima(wbbcf,order=c(0,1,2))
pred.arima = predict(fit.arima,n.ahead=500)
plot(wbbcf,xlim=c(1,800),ylim=range(c(-3.2,1)))
lines(pred.arima$pred,col="red")
lines(pred.arima$pred+1.96*pred.arima$se,col="blue",lty=3)
lines(pred.arima$pred-1.96*pred.arima$se,col="green",lty=3)

Here's the auto.arima() output:

Series: wbbcf 
ARIMA(0,1,2)                    

Coefficients:
          ma1      ma2
      -0.3023  -0.3188
s.e.   0.0397   0.0396

sigma^2 estimated as 0.05788:  log likelihood=3.02
AIC=-0.04   AICc=0.01   BIC=12.98

Is something wrong with what I'm doing with ARIMA?

Answer 1

As @EdM suggests, this may be the best it gets with a sinusoid. Sometimes, a sinusoid just isn't really the underlying relationship that the data is being generated from.

Here are four suggestions that are guaranteed to do at least as good, if not better. I'll assume you've dealt with the underlying time trend and that we're simply discussing the cyclical or seasonal pattern in the data.

First, if you're fitting something like $$\hat{y}_t = c_0 + c_1t + c_2\text{sin}(c_3 2\pi (t + a))$$ to your data instead of something like $$\hat{y}_t = c_0 + c_1t + f(t)\text{sin}(c_3 2\pi (t + a)),$$ you may consider the latter. The function $f(t)$ allows for the amplitude of your sinusoid to vary over its support (i.e. the interval $[\min\{t\}, \max\{t\}]$). This additional flexibility allows for the sinusoid to better accommodate the increased range occurring post-2015. I will acknowledge, however, that sinusoidal modeling can get very complicated for this very reason.

Another solution that would certainly do better than your sine wave in capturing the trend in the data would be to use a spline. I think a cubic spline would be sufficient. You can read more about splines in Chapter 5 of this book. The drawback of this approach is that it does poorly extrapolating into the future periods and it also will require significantly more parameters (which may not be problem judging by how your data looks).

A third and likely more robust solution is to use dummy variables to parametrize whatever "seasons" appear to be apparent in the data. For example, by looking at the 2014 year, we can see that roughly the first four months show an increasing linear trend, the second four months seem flat, and the third four months show a decreasing linear trend. If you include a dummy variable for the slope term of some linear regression of $y$ on $t$ (specifically you would have two of them, perhaps for the second and third trimesters only), you would be able to capture this increase-flat-decrease patter. You would also need to account for the fact that the constant term can different in each season as well. For the linear case, you will at most be estimating six parameters in total.

Finally, we can again notice that the seasonal pattern occurring every four months and difference the data accordingly. That is, instead of working with the series $\{y_t\}_{t=2014}^{t=2015.5}$, you may consider working with $\{y_t - y_{t-\tfrac{1}{3}}\}$, where the $\tfrac{1}{3}$ represents four months or one trimester. The only problem is that you currently don't have enough data to do this and, in general it's difficult to ascertain seasonal or cyclical trends when not enough time has elapsed. Thus, I only recommend this method if you can acquire data going back to at least 2005 or earlier.