library(uGMAR)
context("quantile_residual_tests")

## The results slightly differ depending on whether numerical integration (without the package "gsl")
# or hypergeometric function is used to calculate the quantile residuals. Also, the limited precission
# of floating point numbers (default precission used) might accumulate to differences big enough for
# the tests to fail with too small tolerance (this particularly happens if the parameters contain
# large degree of freedom parameters).

# quantile_residual_tests calls other functions to deal with constraints so they are not considered here
params11t <- c(0.9, 0.92, 1.01, 2.89)
params12 <- c(1.7, 0.85, 0.3, 4.12, 0.73, 1.98, 0.63)
params23 <- c(2.7, 0.8, -0.06, 0.3, 3.5, 0.8, -0.07, 2.6, 7.2, 0.3, -0.01, 0.1, 0.6, 0.25)
params23t <- c(2.54, 0.99, -0.21, 0.36, 5.13, 0.9, -0.24, 1.88, 7.93, 0.2, 0.03, 0.1, 0.53, 0.36, 9, 10, 11)
params12r <- c(2.18, 2.93, 0.81, 0.29, 1.9, 0.64)
params12gs <- c(4.13, 0.73, 1.98, 1.7, 0.85, 0.3, 0.37, 9) # M1=1, M2=1
params13gsr <- c(4.8, 3.31, 3.74, 0.69, 2, 0.19, 0.41, 0.34, 0.3, 9) # M1=2, M2=1

stmar11 <- GSMAR(data=simudata[20:50], p=1, M=1, params=params11t, model="StMAR")
gmar12 <- GSMAR(data=simudata[1:100], p=1, M=2, params=params12, model="GMAR")
gmar23 <- GSMAR(data=simudata[120:200], p=2, M=3, params=params23, model="GMAR")
stmar23 <- GSMAR(data=simudata[120:200], p=2, M=3, params=params23t, model="StMAR")
gmar12r <- GSMAR(data=simudata[100:133], p=1, M=2, params=params12r, model="GMAR", restricted=TRUE)
gstmar12 <- GSMAR(data=simudata[1:50], p=1, M=c(1, 1), params=params12gs, model="G-StMAR")
gstmar13r <- GSMAR(data=simudata[1:50], p=1, M=c(2, 1), params=params13gsr, model="G-StMAR", restricted=TRUE)

set.seed(1); qrt11t <- quantile_residual_tests(stmar11, print_res=FALSE, nsimu=20, lags_ac=c(1, 2, 5, 10))
qrt12 <- quantile_residual_tests(gmar12, print_res=FALSE, nsimu=1, lags_ac=c(1, 2, 5, 10))
qrt23 <- quantile_residual_tests(gmar23, lags_ac=c(1), lags_ch=c(3), print_res=FALSE, nsimu=1)
qrt23t <- quantile_residual_tests(stmar23, lags_ac=c(3), lags_ch=c(2), print_res=FALSE, nsimu=1)
qrt12r <- quantile_residual_tests(gmar12r, lags_ac=c(1, 3), lags_ch=c(1, 3), print_res=FALSE, nsimu=1)
qrt12gs <- quantile_residual_tests(gstmar12, lags_ac=c(2), lags_ch=c(1), print_res=FALSE, nsimu=1)
qrt13gsr <- quantile_residual_tests(gstmar13r, lags_ac=c(1), lags_ch=c(1), print_res=FALSE, nsimu=1)

test_that("quantile residual test for normality works", {
  expect_equal(qrt12gs$norm_res$testStat, 0.5604813, tolerance=1e-3)
  expect_equal(qrt13gsr$norm_res$testStat, 1.663794, tolerance=1e-3)
  expect_equal(qrt11t$norm_res$testStat, 1.082087, tolerance=1e-3)
  expect_equal(qrt12$norm_res$testStat, 1.154846, tolerance=1e-3)
  expect_equal(qrt23$norm_res$testStat, 66.34326, tolerance=1e-3)
  expect_equal(qrt23t$norm_res$testStat, 14.69999, tolerance=1e-3)
  expect_equal(qrt23t$norm_res$pvalue, 0.002091835, tolerance=1e-3)
  expect_equal(qrt23t$norm_res$df, 3, tolerance=1e-4)
  expect_equal(qrt12r$norm_res$testStat, 4.479243, tolerance=1e-3)
  expect_equal(qrt12r$norm_res$pvalue, 0.2141493, tolerance=1e-3)
})

test_that("quantile residuals tests for autocorrelation work", {
  expect_equal(qrt12gs$ac_res$testStat, 0.295045, tolerance=1e-3)
  expect_equal(qrt12gs$ac_res$pvalue, 0.862843, tolerance=1e-3)
  expect_equal(qrt13gsr$ac_res$testStat, 0.009712734, tolerance=1e-3)
  expect_equal(qrt11t$ac_res$testStat, c(0.08001943, 0.15416143, 2.86802409, 5.42951337), tolerance=1e-3)
  expect_equal(qrt11t$ac_res$indStat, c(-0.04784512, 0.05948682, 0.03265780, -0.15573950), tolerance=1e-3)
  expect_equal(qrt12$ac_res$testStat, c(0.2706776, 0.3020141, 1.7811957, 5.8232758), tolerance=1e-3)
  expect_equal(qrt12$ac_res$indStat, c(0.04760758, 0.02580472, 0.04113457, 0.02893367), tolerance=1e-3)
  expect_equal(qrt12$ac_res$pvalue, c(0.6028776, 0.8598416, 0.8785159, 0.8298858), tolerance=1e-3)
  expect_equal(qrt12$ac_res$stdError, c(0.09104275, 0.09432279, 0.08383650, 0.09453681), tolerance=1e-3)
  expect_equal(qrt23$ac_res$testStat, 9.857347, tolerance=1e-3)
  expect_equal(qrt23$ac_res$indStat, 0.2548895, tolerance=1e-3)
  expect_equal(qrt23$ac_res$stdError, 0.08066882, tolerance=1e-3)
  expect_equal(qrt23t$ac_res$testStat, 1.123475, tolerance=1e-3)
  expect_equal(qrt23t$ac_res$pvalue, 0.7714103, tolerance=1e-3)
  expect_equal(qrt23t$ac_res$indStat, -0.01224414, tolerance=1e-3)
  expect_equal(qrt23t$ac_res$stdError, 0.1135299, tolerance=1e-3)
  expect_equal(qrt12r$ac_res$testStat, c(0.4686411, 1.0670959), tolerance=1e-3)
  expect_equal(qrt12r$ac_res$indStat, c(0.1396998, -0.1936938), tolerance=1e-3)
  expect_equal(qrt12r$ac_res$testStat, c(0.4686411, 1.0670959), tolerance=1e-3)
  expect_equal(qrt12$ac_res$stdError, c(0.09104275, 0.09432279, 0.08383650, 0.09453681), tolerance=1e-3)
})

test_that("quantile residual tests for conditional heteroskedasticity work", {
  expect_equal(qrt12gs$ch_res$testStat, 0.8792476, tolerance=1e-3)
  expect_equal(qrt12gs$ch_res$pvalue, 0.3484078, tolerance=1e-3)
  expect_equal(qrt13gsr$ch_res$indStat, -0.2648605, tolerance=1e-3)
  expect_equal(qrt11t$ch_res$testStat, c(0.1806701, 0.9454276, 4.3332202, 221.6827652), tolerance=1e-3)
  expect_equal(qrt11t$ch_res$stdError, c(0.2823567, 0.2931707, 0.2572787, 0.3411960), tolerance=1e-3)
  expect_equal(qrt12$ch_res$testStat, c(0.006418012, 0.056185882, 1.149738575, 7.252740357), tolerance=1e-3)
  expect_equal(qrt12$ch_res$pvalue, c(0.9361478, 0.9722980, 0.9496074, 0.7013887), tolerance=1e-3)
  expect_equal(qrt12$ch_res$indStat, c(0.01729458, -0.04400772, -0.08094851, 0.05611263), tolerance=1e-3)
  expect_equal(qrt12$ch_res$stdError, c(0.2147856, 0.1847421, 0.1136077, 0.1321347), tolerance=1e-3)
  expect_equal(qrt23$ch_res$testStat, 1.525931, tolerance=1e-3)
  expect_equal(qrt23$ch_res$indStat, -0.02697828, tolerance=1e-3)
  expect_equal(qrt23$ch_res$stdError, 0.2194274, tolerance=1e-3)
  expect_equal(qrt23t$ch_res$testStat, 1.51146, tolerance=1e-3)
  expect_equal(qrt23t$ch_res$pvalue, 0.4696677, tolerance=1e-3)
  expect_equal(qrt23t$ch_res$indStat, -0.2311205, tolerance=1e-3)
  expect_equal(qrt23t$ch_res$stdError, 0.2523462, tolerance=1e-3)
  expect_equal(qrt12r$ch_res$testStat, c(0.4177885, 0.3358321), tolerance=1e-3)
  expect_equal(qrt12r$ch_res$testStat, c(0.4177885, 0.3358321), tolerance=1e-3)
  expect_equal(qrt12r$ch_res$indStat, c(-0.18011478, -0.07386617), tolerance=1e-3)
  expect_equal(qrt12r$ch_res$stdError, c(0.2744032, 0.4831642), tolerance=1e-3)
})