beanumber / tidychangepoint

Changepoint detection with a tidy interface
https://beanumber.github.io/tidychangepoint/
Other
0 stars 1 forks source link

discrepancy in computation of log-likelihood for AR(1) model #73

Closed beanumber closed 2 weeks ago

beanumber commented 3 months ago

The computation of the log-likelihood for the trend shift model with AR(1) errors is off by a small amount. However, as in #72 , the calculation is exact for white noise errors. The value of $\sigma^2$ is also correct (to three digits). Even though the difference is small, I'd like to close the gap.

library(tidyverse)
library(tidychangepoint)
library(testthat)
#> 
#> Attaching package: 'testthat'
#> The following object is masked from 'package:dplyr':
#> 
#>     matches
#> The following object is masked from 'package:purrr':
#> 
#>     is_null
#> The following objects are masked from 'package:readr':
#> 
#>     edition_get, local_edition
#> The following object is masked from 'package:tidyr':
#> 
#>     matches

cpts <- c(1700, 1739, 1988)
ids <- time2tau(cpts, substr(time(CET), 1, 4))

trend_ar1 <- fit_lmshift(CET, tau = ids, trends = TRUE, ar1 = TRUE)
expect_equal(round(trend_ar1$sigma_hatsq, 3), 0.290)
expect_equal(round(as.numeric(logLik(trend_ar1)), 2), -288.80)
#> Error: round(as.numeric(logLik(trend_ar1)), 2) not equal to -288.8.
#> 1/1 mismatches
#> [1] -289 - -289 == -0.62
logLik(trend_ar1)
#> 'log Lik.' -289.4169 (df=13)

Created on 2024-04-03 with reprex v2.1.0

beanumber commented 3 months ago

@xuehens do you know what might be going on here?

beanumber commented 2 months ago

I don't know what could be happening here, because the formula is just:

ll <- -(N * log(sigma_hatsq) + N + N * log(2 * pi)) / 2

so the only variable is $\hat{\sigma}^2$, and we already know that is correct to three digits. The value of $\hat{\sigma}^2$ that would produce a log-likelihood of -288.80 is 0.28872, whereas we're getting 0.28971 and the paper reports 0.290.

Interestingly, using $N = 361$ instead of 362 gives the correct value. Is that because you lose an observation by lagging the errors??!

beanumber commented 2 months ago

I think we're within the realm of rounding errors here...