How to calculate the area under each end of a sine curve

Given this data set:

y<-c(-13,16,35,40,28,36,43,33,40,33,22,-5,-27,-31,-29,-25,-26,-31,-26,-24,-25,-29,-23,4)
t<-1:24

My goal is to calculate two areas. The first area would integrate only data from the first part of the curve found above the Zero line. The second area would integrate data from the second part of the curve found below the zero line.

First I would like to fit a sine wave to this data. Using this excellent answer:

https://stats.stackexchange.com/questions/60994/fit-a-sinusoidal-term-to-data

I was able to fit a sine wave (I will be using the periodic with second harmonic which looks to have a better fit)

ssp <- spectrum(y)  
per <- 1/ssp$freq[ssp$spec==max(ssp$spec)]
reslm <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t))
summary(reslm)

rg <- diff(range(y))
plot(y~t,ylim=c(min(y)-0.1*rg,max(y)+0.1*rg))
lines(fitted(reslm)~t,col=4,lty=2)   # dashed blue line is sin fit

# including 2nd harmonic really improves the fit
reslm2 <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t))
summary(reslm2)
lines(fitted(reslm2)~t,col=3)    # solid green line is periodic with second harmonic
abline(h=0,lty=2)

Next I would like to calculate the area under the curve that is only positive, as well as the area under the curve that is exclusively negative. I've had luck looking at similar answers using the AUC functions in the Bolstad2 and Mess packages. But my data points do not fall neatly on zero line, and I do not know how to break up the sine function into areas only above the Zero line and only below the Zero line.

Solution

First things first. To get an exact calculation, you will need to work with the exact function of the 2nd harmonic fourier. Secondly, the beauty of harmonics functions is that they are repetitive. So if you want to find where your function reaches 0, you merely need to expand your interval to so you can be sure to find more than 2 roots.

First we get the exact function from the regression model

fourierfnct <- function(t){
  fnct <- reslm2$coeff[1]+
    reslm2$coeff[2]*sin(2*pi/per*t)+
    reslm2$coeff[3]*cos(2*pi/per*t)+
    reslm2$coeff[4]*sin(4*pi/per*t)+
    reslm2$coeff[5]*cos(4*pi/per*t)
  return(fnct)
}

secondly,you can write a function which can find the roots (where the function is 0). R provides a uniroot function which you can use to find multiple roots in a loop.

manyroots <- function(f,inter,period){
  roots <- array(NA, inter)
  for(i in 1:(length(inter)-1)){
    roots[i] <- tryCatch({
      return_value <- uniroot(f,c(inter[i],inter[i+1]))$root
    }, error = function(err) {
      return_value <- -1
    })
  }
  retroots <- roots[-which(roots==-1)]
  return(retroots)
}

then you simply calculate the roots, and use them to integrate the function across those boundaries.

roots <- manyroots(fourierfnct,seq(0,25),per)
integrate(fourierfnct, roots[1],roots[2])
#300.6378 with absolute error < 3.3e-12
integrate(fourierfnct, roots[2],roots[3])
#-284.6378 with absolute error < 3.2e-12