lubridate,

crayon,

utils,

FKF,

nloptr

nloptr,

BB,

GA,

dfoptim,

pso,

FKF

Suggests:

longmemo,

yardstick,

tidyverse,

BB,

GA,

pso,

dfoptim,

testthat,

knitr,

rmarkdown

export(gg_raw_pgram)

export(ggbr_semipara)

export(version)

importFrom(BB,BBoptim)

importFrom(GA,de)

importFrom(Rsolnp,gosolnp)

importFrom(Rsolnp,solnp)

importFrom(dfoptim,hjkb)

importFrom(dfoptim,nmkb)

importFrom(forecast,forecast)

importFrom(forecast,ggtsdisplay)

importFrom(ggplot2,aes)

importFrom(pracma,deconv)

importFrom(pracma,hessian)

importFrom(pracma,pinv)

importFrom(pso,psoptim)

importFrom(signal,Arma)

importFrom(signal,filter)

importFrom(stats,Box.test)

importFrom(stats,end)

importFrom(stats,frequency)

importFrom(stats,na.pass)

importFrom(stats,optim)

importFrom(stats,optimise)

importFrom(stats,sd)

importFrom(stats,spectrum)

Version 0.9.7 adds the "tsdiag" function for garma models - we just call the base version of tsdiag, so all parameters for that should work.

Version 0.9.7 adds the "tsdiag" function for garma models - this is essentially a copy of the base stats 'tsdiag'

but the fitdf for the Ljung-Box test is set to p+q+k*2.

The override of the ggplot function for garma models has been removed as this was not standard - the correct way to do this is 'autoplot'

so now the autoplot function has the ability to generate forecasts and plot them.

There was a bug with the standard errors which resulted in Nan being returned for some parameters. This has been fixed.

The override of the ggplot function for garma models has been removed as this was not standard - the correct way to do this

is 'autoplot' so now the autoplot function has the ability to generate forecasts and plot them.

A new optimisation method has been added for the "garma" function - this is a genetic algorithm from package GA.

It can be used by specifying opt_method='ga'. For now please treat this as experimental - it has given some interesting

results on some standard problems and these probably need to be checked further before more general use.

Otherwise some redundant code has been removed, and also the restriction that integer differencing be restricted to 1 only.

stop('QML method does not support k>1. It is suggested you try either the CSS or Whittle methods.\n')

if (is.null(opt_method)) {

if (k>=2)opt_method<-'solnp'

elseopt_method<- c('directL','solnp')

opt_method<- c('directL','solnp')

}

for (ominopt_method) {

pars<-numeric(0)

lb<-numeric(0)

ub<-numeric(0)

lb_finite<-numeric(0)

ub_finite<-numeric(0)

mean_methods<- c('QML','CSS')

if (include.mean&method%in%mean_methods) {

mean_y<- mean(y)

pars<- c(pars,mean_y)

lb_finite<- c(lb_finite,ifelse(mean_y<0,2*mean_y,-2*mean_y))

ub_finite<- c(ub_finite,ifelse(mean_y<0,-2*mean_y,2*mean_y))

lb<- c(lb,-Inf)

ub<- c(ub,Inf)

lb<- c(lb,ifelse(mean_y<0,2*mean_y,-2*mean_y))

ub<- c(ub,ifelse(mean_y<0,-2*mean_y,2*mean_y))

}

if (k>0) {# initial parameter estimates for Gegenbauer factors

pars<- c(pars,start_u,start_d)

lb<- c(lb,min_u,ifelse(allow_neg_d,-1,0))

ub<- c(ub,max_u,0.5)

lb_finite<- c(lb_finite,min_u,ifelse(allow_neg_d,-1,0))

ub_finite<- c(ub_finite,max_u,0.5)

}

}

if (p==1&q==0) {# special limits for AR(1)

lb<-c(lb,-1)

if (method%in%methods_to_estimate_var)lb<- c(lb,1e-10)

lb_finite<- c(lb_finite,-1)

if (method%in%methods_to_estimate_var)lb_finite<- c(lb_finite,1e-10)

ub<-c(ub,1)

if (method%in%methods_to_estimate_var)ub<-c(ub,Inf)

ub_finite<- c(ub_finite,1)

if (method%in%methods_to_estimate_var)ub_finite<- c(ub_finite,2*stats::var(y))

if (method%in%methods_to_estimate_var)ub<-c(ub,2*stats::var(y))

}else {

lb<- c(lb,rep(-Inf,p+q))

lb<- c(lb,rep(-10,p+q))

if (method%in%methods_to_estimate_var)lb<- c(lb,1e-10)

ub<- c(ub,rep(Inf,p+q))

if (method%in%methods_to_estimate_var)ub<- c(ub,Inf)

lb_finite<- c(lb_finite,rep(-10,p+q))

if (method%in%methods_to_estimate_var)lb_finite<- c(lb_finite,1e-10)

ub_finite<- c(ub_finite,rep(10,p+q))

if (method%in%methods_to_estimate_var)ub_finite<- c(ub_finite,2*stats::var(y))

ub<- c(ub,rep(10,p+q))

if (method%in%methods_to_estimate_var)ub<- c(ub,2*stats::var(y))

}

}

message<- c()

fcns<-list('CSS'=.css.ggbr.obj,'Whittle'=.whittle.ggbr.obj,'QML'=.qml.ggbr.obj,'WLL'=.wll.ggbr.obj)

fcns<-list('CSS'=.css.ggbr.obj,'Whittle'=.whittle.ggbr.obj,'QML'=.qml.ggbr.obj,'WLL'=.wll.ggbr.obj)

grads<-list('CSS'=NULL,'Whittle'=.whittle.ggbr.grad,'QML'=NULL,'WLL'=NULL)

if (opt_method[[1]]=='best') {

fit<- .best_optim(initial_pars=pars,fcn=fcns[[method]],lb=lb,ub=ub,lb_finite=lb_finite,ub_finite=ub_finite,params=params,control=control)

fit<- .best_optim(initial_pars=pars,fcn=fcns[[method]],grad=grads[[method]],

lb=lb,ub=ub,params=params,control=control)

}else {

fit<- .generic_optim_list(opt_method_list=opt_method,initial_pars=pars,fcn=fcns[[method]],

lb=lb,ub=ub,lb_finite=lb_finite,ub_finite=ub_finite,params=params,control=control)

fit<- .generic_optim_list(opt_method_list=opt_method,initial_pars=pars,fcn=fcns[[method]],grad=grads[[method]],

lb=lb,ub=ub,params=params,control=control)

}

if (fit$convergence==-999) stop('Failed to converge.')

if (fit$convergence==-999) stop('ERROR:Failed to converge.')

if (opt_method[[1]]=='solnp')conv_ok<- (fit$convergence==0)

elseconv_ok<- (fit$convergence>=0)

vcov1<-matrix(nrow=length(fit$par),ncol=length(fit$par))# default

if (conv_ok&method!='WLL'&!is.null(hh)) {

se<-c()# default to set this up in the right environment

h_inv_diag<- diag(inv_hessian<-pracma::pinv(hh))

if (!any(h_inv_diag<0)) {# if hessian doesn't look valid then... generate another one

if (any(h_inv_diag<0)) {# if hessian doesn't look valid then... generate another one

hh<-pracma::hessian(fcns[[method]],fit$par,params=params)

h_inv_diag<- diag(inv_hessian<-pracma::pinv(hh))

}

if (method=='Whittle') {

se<- suppressWarnings(sqrt(2*pi*h_inv_diag))#var(y)/length(y) *

}

if (method=='CSS') {

se<- suppressWarnings(sqrt(h_inv_diag*sigma2*2))

}

if (method=='QML') {

se<- suppressWarnings(sqrt(h_inv_diag*length(y)))

vcov<-inv_hessian*length(y)

vcov1<-inv_hessian*length(y)

}

}

if (method=='WLL') {

se<-rep(NA,length(par))

if (k==1)se[2]<- .wll_d_se(fit$par[1],ss)# result only holds for k=1

vcov<-matrix(nrow=length(fit$par),ncol=length(fit$par))

}

nm<-list()

if (include.mean)nm<- c(nm,'intercept')

coef<- t(matrix(c(temp_coef,temp_se),nrow=n_coef))

colnames(coef)<-nm

rownames(coef)<- c('','s.e.')

colnames(vcov)<- rownames(vcov)<- tail(nm,nrow(vcov))

colnames(vcov1)<- rownames(vcov1)<- tail(nm,nrow(vcov1))

gf<-list()

res<-list('call'= match.call(),

'series'= deparse(match.call()$x),

'coef'=coef,

'var.coef'=vcov,

'var.coef'=vcov1,

'sigma2'=sigma2,

'obj_value'=fit$value,

'loglik'=loglik,

res<- sum(I_f,na.rm=TRUE)#- sum(log(spec_den_inv[spec_den_inv>0]),na.rm=TRUE)

res<- sum(I_f,na.rm=TRUE)- sum(log(spec_den_inv[spec_den_inv>0]),na.rm=TRUE)

return(res)

}

return(res)

}

.whittle.ggbr.grad<-function(par,params) {

n_freq<-length(ss$freq)

freq<-ss$freq[1:n_freq]

spec<-ss$spec[1:n_freq]

u<-fd<- c()

for (k1in seq_len(k)) {

u<- c(u,par[start])

fd<- c(fd,par[start+1])

}

if (p>0)phi_vec<- (-par[start:(start+p-1)])elsephi_vec<-1

if (q>0)theta_vec<- (par[(p+start):length(par)])elsetheta_vec<-1

cos_2_pi_f<- cos(2.0*pi*freq)

if (p>0)mod_phi<- .a_fcn(phi_vec,freq)

if (q>0)mod_theta<- .a_fcn(theta_vec,freq)

for (k1in seq_len(k)) {

u_k<-u[k1]

fd_k<-fd[k1]

temp<- (4.0*((cos_2_pi_f-u_k)^2.0))

spec_den_inv<-spec_den_inv* (temp^fd_k)

}

spec_den_inv[spec_den_inv==0]<-NA

grad<- c()

for (k1in seq_len(k)) {

grad_u<- sum(2*fd[k1]/(cos_2_pi_f-u[k1])*(1.0-I_f),na.rm=TRUE)

grad_fd<- sum(-log((cos_2_pi_f-u[k1])*(cos_2_pi_f-u[k1])*4)*(1.0-I_f),na.rm=TRUE)

grad<- c(grad,grad_u,grad_fd)

}

if (p>0) {

phi_cos_vec<- .a_fcn_cos(phi_vec,freq)

phi_sin_vec<- .a_fcn_sin(phi_vec,freq)

for (p1in seq_len(p)) {

grad_phi<-2.0*sum((phi_cos_vec*cos(p1*2*pi*freq)+phi_sin_vec*sin(p1*2*pi*freq))/mod_phi* (1.0-I_f),na.rm=TRUE)

grad<- c(grad,grad_phi)

}

}

if (q>0) {

theta_cos_vec<- .a_fcn_cos(theta_vec,freq)

theta_sin_vec<- .a_fcn_sin(theta_vec,freq)

for (q1in seq_len(q)) {

grad_theta<-2.0*sum((theta_cos_vec*cos(q1*2*pi*freq)+theta_sin_vec*sin(q1*2*pi*freq))/mod_theta* (I_f+1.0),na.rm=TRUE)

grad<- c(grad,grad_theta)

}

}

return(grad)

}