functions.R

lag.plot1=function(series,max.lag=1,corr=TRUE,smooth=FALSE){ 
   name1=paste(deparse(substitute(series)),"(t-",sep="")
   name2=paste(deparse(substitute(series)),"(t)",sep="")
   data1=as.ts(series)
   max.lag=as.integer(max.lag)
   prow=ceiling(sqrt(max.lag))
   pcol=ceiling(max.lag/prow)
   a=acf(series,max.lag,plot=FALSE)$acf[-1]
   old.par <- par(no.readonly = TRUE)
   par(mfrow=c(prow,pcol), mar=c(2.5, 4, 2.5, 1), cex.main=1.1, font.main=1)
  for(h in 1:max.lag){                       
   plot(lag(series,-h), data1, xy.labels=FALSE, main=paste(name1,h,")",sep=""), ylab=name2, xlab="") 
    if (smooth==TRUE) 
    lines(lowess(ts.intersect(lag(series,-h),series)[,1],
                 ts.intersect(lag(series,-h),series)[,2]), col="red")
    if (corr==TRUE)
    legend("topright", legend=round(a[h], digits=2), text.col ="blue", bg="white", x.intersp=0)
   on.exit(par(old.par))
   }
}

lag.plot2=function(series1,series2,max.lag=0,corr=TRUE,smooth=FALSE){ 
   name1=paste(deparse(substitute(series1)),"(t-",sep="")
   name2=paste(deparse(substitute(series2)),"(t)",sep="")
   data1=as.ts(series1)
   data2=as.ts(series2)
   max.lag=as.integer(max.lag)
   m1=max.lag+1
   prow=ceiling(sqrt(m1))
   pcol=ceiling(m1/prow)
   a=ccf(series1,series2,max.lag,plot=FALSE)$acf
   old.par <- par(no.readonly = TRUE)
   par(mfrow=c(prow,pcol), mar=c(2.5, 4, 2.5, 1), cex.main=1.1, font.main=1)
   for(h in 0:max.lag){                   
   plot(lag(series1,-h), series2, xy.labels=FALSE, main=paste(name1,h,")",sep=""), ylab=name2, xlab="") 
    if (smooth==TRUE) 
    lines(lowess(ts.intersect(lag(series1,-h),series2)[,1],
                 ts.intersect(lag(series1,-h),series2)[,2]), col="red")
    if (corr==TRUE)
    legend("topright", legend=round(a[m1-h], digits=2), text.col ="blue", bg="white", x.intersp=0)             
   on.exit(par(old.par))
   }
}

acf2=function(series,max.lag=NULL){
  num=length(series)
  if (is.null(max.lag)) max.lag=ceiling(10+sqrt(num))
  if (max.lag > (num-1)) stop("Number of lags exceeds number of observations")
  ACF=acf(series, max.lag, plot=FALSE)$acf[-1]
  PACF=pacf(series, max.lag, plot=FALSE)$acf
  LAG=1:max.lag/frequency(series)
  minA=min(ACF)
  minP=min(PACF)
  U=2/sqrt(num)
  L=-U
  minu=min(minA,minP,L)-.01
  old.par <- par(no.readonly = TRUE)
  par(mfrow=c(2,1), mar = c(3,3,2,0.8),
    oma = c(1,1.2,1,1), mgp = c(1.5,0.6,0))
  plot(LAG, ACF, type="h",ylim=c(minu,1), 
    main=paste("Series: ",deparse(substitute(series))))
    abline(h=c(0,L,U), lty=c(1,2,2), col=c(1,4,4))
  plot(LAG, PACF, type="h",ylim=c(minu,1))
    abline(h=c(0,L,U), lty=c(1,2,2), col=c(1,4,4))
  on.exit(par(old.par))  
  ACF<-round(ACF,2); PACF<-round(PACF,2)    
  return(cbind(ACF, PACF)) 
  }

sarima = function(xdata,p,d,q,P=0,D=0,Q=0,S=-1,details=TRUE,tol=sqrt(.Machine$double.eps),no.constant=FALSE)
{ 
  n = length(xdata)
  constant = 1:n   
  xmean = rep(1,n) 
  trc = ifelse(details==TRUE, 1, 0)
  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, 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,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), 
              optim.control=list(trace=trc,REPORT=1,reltol=tol))
#
#  replace tsdiag with a better version
    old.par <- par(no.readonly = TRUE)
    layout(matrix(c(1,2,4, 1,3,4), nc=2))
    rs <- fitit$residuals
    stdres <- rs/sqrt(fitit$sigma2)
    num <- sum(!is.na(rs))
     plot.ts(stdres,  main = "Standardized Residuals", ylab = "")
    alag <- 10+sqrt(num)
    ACF = acf(rs, alag, plot=FALSE, na.action = na.pass)$acf[-1] 
    LAG = 1:alag/frequency(xdata)
    L=2/sqrt(num)
     plot(LAG, ACF, type="h", ylim=c(min(ACF)-.1,min(1,max(ACF+.4))), main = "ACF of Residuals")
     abline(h=c(0,-L,L), lty=c(1,2,2), col=c(1,4,4))  
    qqnorm(stdres, main="Normal Q-Q Plot of Std Residuals"); qqline(stdres, col=4) 
    nlag <- ifelse(S<4, 20, 3*S)
    ppq <- p+q+P+Q
    pval <- numeric(nlag)
    for (i in (ppq+1):nlag) {u <- Box.test(rs, i, type = "Ljung-Box")$statistic
                             pval[i] <- pchisq(u, i-ppq, lower=FALSE)}            
     plot( (ppq+1):nlag, pval[(ppq+1):nlag], xlab = "lag", ylab = "p value", ylim = c(0, 
        1), main = "p values for Ljung-Box statistic")
     abline(h = 0.05, lty = 2, col = "blue")  
    on.exit(par(old.par))    
#  end new tsdiag
#
  k = length(fitit$coef)
  BIC = log(fitit$sigma2)+(k*log(n)/n)
  AICc = log(fitit$sigma2)+((n+k)/(n-k-2))
  AIC = log(fitit$sigma2)+((n+2*k)/n)
  innov<<-fitit$resid
  list(fit=fitit, AIC=AIC, AICc=AICc, BIC=BIC)
}

sarima.for=function(xdata,n.ahead,p,d,q,P=0,D=0,Q=0,S=-1,tol=sqrt(.Machine$double.eps),no.constant=FALSE){ 
  xname=deparse(substitute(xdata))
  xdata=as.ts(xdata) 
  n=length(xdata)
  constant=1:n
  xmean=rep(1,n)
  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, optim.control=list(reltol=tol));
    nureg=matrix(1,n.ahead,1)        
} 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,optim.control=list(reltol=tol));
    nureg=(n+1):(n+n.ahead)       
} else { fitit=arima(xdata, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S), 
            optim.control=list(reltol=tol));
          nureg=NULL   
}     
#--
 fore=predict(fitit, n.ahead, newxreg=nureg)  
#-- graph:
  U = fore$pred + 2*fore$se
  L = fore$pred - 2*fore$se
   a=max(1,n-100)
  minx=min(xdata[a:n],L)
  maxx=max(xdata[a:n],U)
   t1=xy.coords(xdata, y = NULL)$x 
   if(length(t1)<101) strt=t1[1] else strt=t1[length(t1)-100]
   t2=xy.coords(fore$pred, y = NULL)$x 
   endd=t2[length(t2)]
   xllim=c(strt,endd)
  ts.plot(xdata,fore$pred,col=1:2, xlim=xllim, ylim=c(minx,maxx), ylab=xname) 
  lines(fore$pred, col="red", type="p")
  lines(U, col="blue", lty="dashed")
  lines(L, col="blue", lty="dashed") 
#
  return(fore)
}


spec.arma <- function(ar=0,ma=0,var.noise=1,n.freq=500, ...)
{ 
    # check causality
     ar.poly <- c(1, -ar)
     z.ar <- polyroot(ar.poly)
     if(any(abs(z.ar) <= 1)) cat("WARNING: Model Not Causal", "\n")  
    # check invertibility
     ma.poly <- c(1, ma)
     z.ma <- polyroot(ma.poly)
     if(any(abs(z.ma) <= 1)) cat("WARNING: Model Not Invertible", "\n")
     if(any(abs(z.ma) <= 1) || any(abs(z.ar) <= 1) ) stop("Try Again")
    #
    ar.order <- length(ar)
    ma.order <- length(ma) 
    # check (near) parameter redundancy [i.e. are any roots (approximately) equal]  
       for (i in 1:ar.order) {
       if ( (ar == 0 & ar.order == 1) || (ma == 0 & ma.order ==1) ) break
       if(any(abs(z.ar[i]-z.ma[1:ma.order]) < 1e-03)) {cat("WARNING: Parameter Redundancy", "\n"); break}
       }
    #
    freq <- seq.int(0, 0.5, length.out = n.freq)
            cs.ar <- outer(freq, 1:ar.order, function(x, y) cos(2 * 
                pi * x * y)) %*% ar
            sn.ar <- outer(freq, 1:ar.order, function(x, y) sin(2 * 
                pi * x * y)) %*% ar
            cs.ma <- outer(freq, 1:ma.order, function(x, y) cos(2 * 
                pi * x * y)) %*% -ma
            sn.ma <- outer(freq, 1:ma.order, function(x, y) sin(2 * 
                pi * x * y)) %*% -ma                      
    spec <- var.noise*((1 - cs.ma)^2 + sn.ma^2)/((1 - cs.ar)^2 + sn.ar^2)
    spg.out <- list(freq=freq, spec=spec)
    class(spg.out) <- "spec"
    plot(spg.out, ci=0, main="", ...)
    return(invisible(spg.out))
}

find.best.arima <- function(Data, maxord=c(2,1,1,2,1,1))
{
	if(!is.ts(Data))
		stop

	Data <- na.trim(Data)

	if(length(Data) < 5 * frequency(Data))
		stop

	best.aic <- 1e8
	n <- length(Data)
	for(p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
	{
		for(P in 0:maxord[4]) for (D in 0:maxord[5]) for(Q in 0:maxord[6])
		{
			fit <- arima(Data, order=c(p,q,d),seas=list(order=c(P,D,Q),
				frequency(Data)), method='CSS')
			fit.aic <- -2*fit$loglik + (log(n) + 1) * length(fit$coef)
			if(fit.aic < best.aic)
			{
				best.aic <- fit.aic
				best.fit <- fit
				p <- p
				d <- d
				q <- q
				P <- P
				D <- D
				Q <- Q
			}
		}
	}	
	return(list(p,d,q,P,D,Q))
}

cat("  itall has been installed", "\n")