Movatterモバイル変換


[0]ホーム

URL:


Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

Sign up
Appearance settings

Commit07ec33c

Browse files
committed
function name changes
1 parent0ce04d7 commit07ec33c

File tree

7 files changed

+64
-36
lines changed

7 files changed

+64
-36
lines changed

‎NAMESPACE‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,11 +15,13 @@ S3method(residuals,garma_model)
1515
S3method(summary,garma_model)
1616
S3method(tsdiag,garma_model)
1717
S3method(vcov,garma_model)
18+
export(calc_garma_spec_den_inv)
1819
export(extract_arma)
1920
export(garma)
2021
export(garma_ggtsdisplay)
2122
export(gg_raw_pgram)
2223
export(ggbr_semipara)
24+
export(gof)
2325
export(version)
2426
importFrom(BB,BBoptim)
2527
importFrom(GA,de)

‎R/garma-package.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212
#' @importFrom pso psoptim
1313
#' @importFrom dfoptim hjkb nmkb
1414
#' @importFrom GA de
15-
#' @importFrom pracma pinv hessian conv deconv grad psi
15+
#' @importFrom pracma pinv hessian conv deconv grad psi digamma
1616
#' @importFrom signal Arma filter
1717
#' @importFrom zoo zoo index
1818
#' @importFrom tswge factor.wge

‎R/garma_main.R‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -269,7 +269,7 @@ garma<-function(x,
269269
message<- c()
270270

271271
# Optimisation functions for each method
272-
fcns<-list('CSS'=.css.ggbr.obj,'Whittle'=.whittle.ggbr.obj.short,'QML'=.qml.ggbr.obj,'WLL'=.wll.ggbr.obj)
272+
fcns<-list('CSS'=.css.ggbr.obj,'Whittle'=.whittle.garma.obj.short,'QML'=.qml.ggbr.obj,'WLL'=.wll.ggbr.obj)
273273
if (opt_method[[1]]=='best') {
274274
fit<- .best_optim(initial_pars=pars,fcn=fcns[[method]],
275275
lb=lb,ub=ub,params=params,control=control)
@@ -288,7 +288,7 @@ garma<-function(x,
288288
}
289289
elseif (method=='QML')sigma2<- .qml.ggbr.se2(fit$par,params=params)
290290
elseif (method=='CSS')sigma2<-fit$value[length(fit$value)]/length(y)# Chung (1996)
291-
elseif (method=='Whittle')sigma2<- .whittle.ggbr.obj.short(fit$par,params)# GHR 2001.
291+
elseif (method=='Whittle')sigma2<- .whittle.garma.obj.short(fit$par,params)# GHR 2001.
292292

293293
# log lik
294294
loglik<-numeric(0)
@@ -297,7 +297,7 @@ garma<-function(x,
297297
if (method=='QML')
298298
loglik<--fit$value[length(fit$value)]
299299
if (method=='Whittle')
300-
loglik<--0.5*(2*length(y)*log(2*pi)+ .whittle.ggbr.obj(fit$par,params))#refer GHR 2001, Whittle (1953) thm 6.
300+
loglik<--0.5*(2*length(y)*log(2*pi)+ .whittle.garma.obj(fit$par,params))#refer GHR 2001, Whittle (1953) thm 6.
301301

302302
se<-numeric(length(fit$par))
303303

@@ -315,7 +315,7 @@ garma<-function(x,
315315
h_inv_diag<- diag(inv_hessian<-pracma::pinv(hh))
316316
if (method=='Whittle') {
317317
# Next line from GHR (2001) thm 3.2
318-
omega<- .whittle_omega(fit$par,params)
318+
omega<- .whittle_garma_omega(fit$par,params)
319319
vcov1<-pracma::pinv(omega)/length(y)
320320
se<- suppressWarnings(sqrt(diag(vcov1)))
321321
}
@@ -345,7 +345,7 @@ garma<-function(x,
345345
temp_se<-se[1:n_coef]
346346
if (include.mean&!method%in%mean_methods) {# then add an Intercept anyway; use Kunsch 1987 result for se
347347
temp_coef<- c(mean(y),temp_coef)
348-
mean_se<-sigma2* .spec_den_0(fit$par,params)/(2*pi)# Chung (1996)
348+
mean_se<-sigma2* .garma_spec_den_0(fit$par,params)/(2*pi)# Chung (1996)
349349
temp_se<- c(mean_se,temp_se)
350350
n_coef<-n_coef+1
351351
}

‎R/garma_whittle.R‎

Lines changed: 12 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -4,13 +4,13 @@
44
# @param par - the parameters to evaluate the function at
55
# @param params - other parameters - including the p, q, k, and scale parameters and (ss) the spectrum .
66
# @return The value of the objective at the point par.
7-
.whittle.ggbr.obj<-function(par,params) {
7+
.whittle.garma.obj<-function(par,params) {
88
# objective function to be minimised for Whittle estimates
99
ss<-params$ss
1010
spec<-ss$spec
1111
n_freq<- length(spec)
1212

13-
spec_den_inv<- .calc_spec_den_inv(par,params)
13+
spec_den_inv<- .calc_garma_spec_den_inv(par,params)
1414

1515
## I have checked the source for for "spectrum" which is used to calculate the "spec"
1616
## The I/f calc includes spectrum and the "f" we calc as above.
@@ -23,14 +23,14 @@
2323
return(res)
2424
}
2525

26-
.whittle.ggbr.obj.short<-function(par,params) {
26+
.whittle.garma.obj.short<-function(par,params) {
2727
# Just the first part of the objective function.
2828
# This is the "S" function of Giraitis, Hidalgo & Robinson 2001.
2929
ss<-params$ss
3030
spec<-ss$spec
3131
n_freq<- length(spec)
3232

33-
spec_den_inv<- .calc_spec_den_inv(par,params)
33+
spec_den_inv<- .calc_garma_spec_den_inv(par,params)
3434

3535
## I have checked the source for for "spectrum" which is used to calculate the "spec"
3636
## The I/f calc includes spectrum and the "f" we calc. f includes the factor 1/(2 pi)
@@ -42,7 +42,10 @@
4242
return(res)
4343
}
4444

45-
.calc_spec_den_inv<-function(par,params) {
45+
## @export
46+
## calc_garma_spec_den_inv<-function(par,params) return(.calc_garma_spec_den_inv(par,params))
47+
48+
.calc_garma_spec_den_inv<-function(par,params) {
4649
# This is the "k" function of Giraitis, Hidalgo & Robinson 2001.
4750
# true spec den is this times sigma^2/(2*pi)
4851
ss<-params$ss
@@ -85,7 +88,7 @@
8588
}
8689

8790

88-
.spec_den_0<-function(par,params) {
91+
.garma_spec_den_0<-function(par,params) {
8992
# Just the first part of the objective function - used for calculating the variance of the process.
9093
ss<-params$ss
9194
p<-params$p
@@ -119,16 +122,16 @@
119122
}
120123

121124
# for vcov calcs...
122-
.whittle_omega<-function(par,params) {
125+
.whittle_garma_omega<-function(par,params) {
123126
# This returns the estimate of the "omega" matrix from GHR (2001). Used for determining se's.
124127

125128
calc_grad_log_k<-function(par,params) {
126129
calc_log_k<-function(par,params) {
127-
k<- .calc_spec_den_inv(par,params)
130+
k<- .calc_garma_spec_den_inv(par,params)
128131
return(log(k))
129132
}
130133
calc_log_k_j<-function(par,params,j) {
131-
k<- .calc_spec_den_inv(par,params)
134+
k<- .calc_garma_spec_den_inv(par,params)
132135
return(log(k[j]))
133136
}
134137

‎R/tsdiag.R‎

Lines changed: 21 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,6 @@
44
#' This function is copied from stats::tsdiag but modifies the fit_df for the Ljung-Box test for use with garma models.
55
#'
66
#' @param object (garma_model) The garma_model to produce the diagnostic plots for.
7-
#' @param gof.lag (int) the maximum number of lags for a Portmanteau goodness-of-fit test.
87
#' @param ... further arguments to be passed to particular methods.
98
#' @return None. Diagnostics are generated.
109
#' @seealso The stats package tsdiag function: \url{https://stat.ethz.ch/R-manual/R-patched/library/stats/html/tsdiag.html}.
@@ -29,21 +28,24 @@ tsdiag.garma_model<-function(object, ...) {
2928
main= paste(object$series," - ACF of Residuals"),
3029
sub=titles$sub,
3130
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')
4531

46-
# Fay & Philippe (2002) spectral test for gof.
32+
gof(object)
33+
}
34+
35+
#' Goodness-of-Fit test for a garma_model.
36+
#'
37+
#' Provides a goodness-of-fit test for a GARMA Model, using the method of Fay & Philippe (2002).
38+
#'
39+
#' @param object (garma_model) The garma_model to produce the diagnostic plots for.
40+
#' @return None.
41+
#' @examples
42+
#' data(AirPassengers)
43+
#' ap <- as.numeric(diff(AirPassengers,12))
44+
#' mdl <- garma(ap,order=c(9,1,0),k=0,method='CSS',include.mean=FALSE)
45+
#' gof(mdl)
46+
#' @export
47+
gof<-function(object) {
48+
# Fay & Philippe (2002) spectral test for gof==goodness-of-fit.
4749
ss<-object$spectrum
4850
n_freq<- length(ss$spec)
4951
mean_methods<- c('QML','CSS')
@@ -53,13 +55,13 @@ tsdiag.garma_model<-function(object, ...) {
5355
est_mean=ifelse(object$method%in%mean_methods,TRUE,FALSE),
5456
scale=sd(object$y),m_trunc=object$m_trunc)
5557

56-
spec_den_inv<- .calc_spec_den_inv(object$obj_par,params)
58+
spec_den_inv<- .calc_garma_spec_den_inv(object$obj_par,params)
5759

5860
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+
s<- log(mean(I_f,na.rm=TRUE))- mean(log(I_f),na.rm=TRUE)+pracma::digamma(1)# (-digamma(1)) == Euler-Mascheroni constant
62+
s_sigma2<-pracma::psi(1,1)/sqrt(n_freq)
6163
s2<-s/s_sigma2
6264
p<-1-pnorm(s2)
6365

64-
cat(sprintf('Fay-Philippe Goodness of Fit test.\n\nz-Statistic: %0.4f\np-value: %0.4f\n',s2,p))
66+
cat(sprintf('Fay-Philippe Goodness of Fit test.\n\nz-statistic: %0.4f\np-value: %0.4f\n',s2,p))
6567
}

‎man/gof.Rd‎

Lines changed: 23 additions & 0 deletions
Some generated files are not rendered by default. Learn more aboutcustomizing how changed files appear on GitHub.

‎man/tsdiag.garma_model.Rd‎

Lines changed: 0 additions & 2 deletions
Some generated files are not rendered by default. Learn more aboutcustomizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)

[8]ページ先頭

©2009-2025 Movatter.jp