@@ -19,25 +19,47 @@ tsdiag.garma_model<-function(object, gof.lag = 10, ...) {
1919titles <- .generate_default_plot_title(object ,h = 0 )
2020
2121# # plot standardized residuals, acf of residuals, Ljung-Box p-values
22- oldpar <- par(mfrow = c(3 ,1 ))
22+ oldpar <- par(mfrow = c(2 ,1 ))
2323 on.exit(par(oldpar ))
2424rs <- object $ residuals
2525stdres <- rs / sqrt(object $ sigma2 )
26- plot(stdres ,type = " h" main = " Standardized Residuals" ylab = " "
26+ plot(stdres ,type = " h" main = paste( object $ series , " - Standardized Residuals" ) ,ylab = " "
2727 abline(h = 0 )
28- acf(as.numeric(object $ residuals ),plot = TRUE ,main = " ACF of Residuals" na.action = na.pass )
29- # pacf(as.numeric(object$residuals), plot = TRUE, main = "PACF of Residuals", na.action = na.pass)
30- nlag <- gof.lag
31- pval <- rep(NA ,nlag )
32- n_param <- object $ order [1 ]+ object $ order [3 ]+ object $ k
33- for (i in 1L : nlag )pval [i ]<- Box.test(rs ,i ,type = " Ljung-Box" fitdf = ifelse(i > n_param ,n_param ,0 ))$ p.value
34- plot(1L : nlag ,pval ,xlab = " Lag" ylab = " p value" ylim = c(0 ,1 ),xaxp = c(1 ,nlag ,nlag - 1 ),
35- main = " p values for Ljung-Box statistic" sub = titles $ sub )
36- abline(h = 0.05 ,lty = 2 ,col = " blue"
28+ acf(as.numeric(object $ residuals ),plot = TRUE ,
29+ main = paste(object $ series ," - ACF of Residuals"
30+ sub = titles $ sub ,
31+ na.action = na.pass )
32+ # nlag <- gof.lag
33+ # pval <- rep(NA,nlag)
34+ # n_param <- object$order[1]+object$order[3]+object$k
35+ # for(i in 1L:nlag) pval[i] <- Box.test(rs, i, type="Ljung-Box", fitdf=ifelse(i>n_param,n_param,0))$p.value
36+ # plot(1L:nlag, pval, xlab = "Lag", ylab = "p value", ylim = c(0,1), xaxp=c(1,nlag,nlag-1),
37+ # main = "p values for Ljung-Box statistic", sub=titles$sub)
38+ # abline(h = 0.05, lty = 2, col = "blue")
39+ #
40+ # cat('NOTE: The degrees of freedom for the Ljung-Box statistic are determined differently\n')
41+ # cat(' to the standard R "tsdiag". An adjustment is made for the number of model\n')
42+ # cat(' parameters as per Ljung & Box (1978).\n')
43+ # cat(' Note also that the Ljung-Box test has not been formally established as valid\n')
44+ # cat(' for GARMA processes.\n')
3745
38- cat(' NOTE: The degrees of freedom for the Ljung-Box statistic are determined differently\n '
39- cat(' to the standard R "tsdiag". An adjustment is made for the number of model\n '
40- cat(' parameters as per Ljung & Box (1978).\n '
41- cat(' Note also that the Ljung-Box test has not been formally established as valid\n '
42- cat(' for GARMA processes.\n '
46+ # Fay & Philippe (2002) spectral test for gof.
47+ ss <- object $ spectrum
48+ n_freq <- length(ss $ spec )
49+ mean_methods <- c(' QML' ' CSS'
50+ params <- list (y = object $ y ,orig_y = object $ diff_y ,ss = ss ,p = object $ order [1 ],q = object $ order [3 ],d = object $ order [2 ],k = object $ k ,
51+ include.mean = object $ include.mean ,include.drift = object $ include.drift ,drift = object $ model $ drift ,
52+ method = object $ method ,
53+ est_mean = ifelse(object $ method %in% mean_methods ,TRUE ,FALSE ),
54+ scale = sd(object $ y ),m_trunc = object $ m_trunc )
55+
56+ spec_den_inv <- .calc_spec_den_inv(object $ obj_par ,params )
57+
58+ I_f <- ss $ spec * spec_den_inv / (2 * pi * n_freq )
59+ s <- log(mean(I_f ,na.rm = TRUE ))- mean(log(I_f ),na.rm = TRUE )- digamma(1 )# (-digamma(1)) == Euler-Mascheroni constant
60+ s_sigma2 <- pracma :: psi(1 ,1 )
61+ s2 <- s / s_sigma2
62+ p <- 1 - pnorm(s2 )
63+
64+ cat(sprintf(' Fay-Philippe Goodness of Fit test.\n\n z-Statistic: %0.4f\n p-value: %0.4f\n ' s2 ,p ))
4365}