library(ggplot2)

x <- reefreviews[complete.cases(reefreviews), ]

# Line plot review ratings
ggplot(data=x, aes(x=month, y=rating)) +
    geom_point()

#autocorrelation test

acf(x$rating,lag.max = 60)
pacf(x$rating,lag.max = 60)

acf(x$reviews,lag.max = 60)
pacf(x$reviews,lag.max = 60)

#de-seasonalize rating time series with bfast

library(bfast)

myts <- ts(x$rating, start=c(2008, 1), end=c(2020, 12), frequency=12)

plot(myts, ylab="NDVI")

(rdist <- 10/length(myts))

reefratings <- bfast(myts,h=rdist, season="harmonic", max.iter=1,breaks=2)

plot(reefratings, main="")

#de-seasonalize number of reviews time series with bfast

myts2 <- ts(x$reviews, start=c(2008, 1), end=c(2020, 12), frequency=12)

plot(myts2, ylab="NDVI")

(rdist <- 10/length(myts2))

reeffrequency <- bfast(myts2,h=rdist, season="harmonic", max.iter=1,breaks=2)

plot(reeffrequency, main="")


#de-seasonalize rating time series with bfast

y <- forestreviews[complete.cases(forestreviews), ]

library(bfast)

myts3 <- ts(y$rating, start=c(2008, 1), end=c(2020, 12), frequency=12)

plot(myts3, ylab="NDVI")

(rdist <- 10/length(myts3))

forestratings <- bfast(myts3,h=rdist, season="harmonic", max.iter=1,breaks=2)

plot(forestratings,main="")

#de-seasonalize number of reviews time series with bfast

myts4 <- ts(y$reviews, start=c(2008, 1), end=c(2020, 12), frequency=12)

plot(myts4, ylab="NDVI")

(rdist <- 10/length(myts4))

forestfrequency <- bfast(myts4,h=rdist, season="harmonic", max.iter=1,breaks=2)

plot(forestfrequency, main="")