Nothing
sarima <-
function(xdata,p,d,q,P=0,D=0,Q=0,S=-1,details=TRUE,xreg=NULL,Model=TRUE,
fixed=NULL,tol=sqrt(.Machine$double.eps),no.constant=FALSE, ...)
{
trans = ifelse (is.null(fixed), TRUE, FALSE)
trc = ifelse(details, 1, 0)
n = length(xdata)
if (is.null(xreg)) {
constant = 1:n
xmean = rep(1,n)
if (no.constant) xmean=NULL
if (d==0 & D==0) {
fitit = arima(xdata, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S),
xreg=xmean, include.mean=FALSE, fixed=fixed, transform.pars=trans,
optim.control=list(trace=trc, REPORT=1, reltol=tol))
} else if (xor(d==1, D==1) & no.constant==FALSE) {
fitit = arima(xdata, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S),
xreg=constant, fixed=fixed, transform.pars=trans,
optim.control=list(trace=trc, REPORT=1, reltol=tol))
} else fitit = arima(xdata, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S),
include.mean=!no.constant, fixed=fixed, transform.pars=trans,
optim.control=list(trace=trc, REPORT=1, reltol=tol))
}
if (!is.null(xreg)) {
fitit = arima(xdata, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S),
xreg=xreg, fixed=fixed, transform.pars=trans,
optim.control=list(trace=trc, REPORT=1, reltol=tol))
}
if (is.null(fixed)) {coefs=fitit$coef
} else { coefs = fitit$coef[is.na(fixed)]
}
k = length(coefs)
n = fitit$nobs # effective sample size
dfree = n-k
t.value = coefs/sqrt(diag(fitit$var.coef))
p.two = stats::pf(t.value^2, df1=1, df2=dfree, lower.tail = FALSE)
ttable = cbind(Estimate=coefs, SE=sqrt(diag(fitit$var.coef)), t.value, p.value=p.two)
ttable = round(ttable,4)
BIC = stats::BIC(fitit)/n
AIC = stats::AIC(fitit)/n
AICc = (n*AIC + ( (2*k^2+2*k)/(n-k-1) ))/n
# print results
cat('<><><><><><><><><><><><><><>')
cat('\n','\n')
cat('Coefficients:', '\n')
print(ttable)
cat('\n')
cat('sigma^2 estimated as', fitit$sigma2, 'on', dfree, 'degrees of freedom', '\n','\n')
cat('AIC =', AIC, ' AICc =', AICc, ' BIC =', BIC, '\n', '\n')
out = list(fit=fitit, degrees_of_freedom=dfree, ttable=ttable,
ICs=c(AIC=AIC, AICc=AICc, BIC=BIC))
############################
# replace tsdiag with a better version
if(details){
old.par <- par(no.readonly = TRUE)
layout(matrix(c(1,2,4, 1,3,4), ncol=2))
rs <- fitit$residuals
stdres <- rs/sqrt(fitit$sigma2)
num <- sum(!is.na(rs))
tsplot(stdres, main = "Standardized Residuals", ylab = "", ...)
if(Model){
if (S<0) {title(bquote('Model: ('~.(p)*','*.(d)*','*.(q)~')'), adj=0) }
else {title(bquote('Model: ('~.(p)*','*.(d)*','*.(q)~')'~'\u00D7'~'('~.(P)*','*.(D)*','*.(Q)~')'[~.(S)]),
adj=0) }
}
alag <- max(10+sqrt(num), 3*S)
acf1(rs, alag, main = "ACF of Residuals", ...)
u = qqnorm(stdres, plot.it=FALSE)
lwr = min(-4, min(stdres, na.rm=TRUE)); upr = max(4, max(stdres), na.rm=TRUE)
tsplot(u$x, u$y, type='p', ylim=c(lwr,upr), ylab="Sample Quantiles", xlab="Theoretical Quantiles",
main="Normal Q-Q Plot of Std Residuals", ...)
################ qq error bnds ###########
sR <- !is.na(stdres)
ord <- order(stdres[sR])
ord.stdres <- stdres[sR][ord]
PP <- stats::ppoints(num)
z <- stats::qnorm(PP)
y <- stats::quantile(ord.stdres, c(.25,.75), names = FALSE, type = 7, na.rm = TRUE)
x <- stats::qnorm(c(.25,.75))
b <- diff(y)/diff(x)
a <- y[1L] - b * x[1L]
abline(a,b,col=4) #qqline
SE <- (b/dnorm(z))*sqrt(PP*(1-PP)/num)
qqfit <- a + b*z
U <- qqfit+3.9*SE # puts .00005 in tails
L <- qqfit-3.9*SE
z[1]=z[1]-.1 # extend plot -- misses the end otherwise
z[length(z)]= z[length(z)]+.1
xx <- c(z, rev(z))
yy <- c(L, rev(U))
polygon(xx, yy, border=NA, col=gray(.5, alpha=.2) )
############ end qq error bnds ##########################
nlag <- ifelse(S<7, 20, 3*S)
ppq <- p+q+P+Q - sum(!is.na(fixed)) # decrease by number of fixed parameters
if (nlag < ppq + 8) {nlag = ppq + 8}
pval <- numeric(nlag)
for (i in (ppq+1):nlag) {u <- stats::Box.test(rs, i, type = "Ljung-Box")$statistic
pval[i] <- stats::pchisq(u, i-ppq, lower.tail=FALSE)}
tsplot( (ppq+1):nlag, pval[(ppq+1):nlag], type='p', xlab = "LAG (H)", ylab = "p value",
ylim = c(-.14, 1), main = "p values for Ljung-Box statistic", ...)
abline(h = 0.05, lty = 2, col = 4)
on.exit(par(old.par))
} # end new tsdiag
invisible(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.