sarima=function(data,p,d,q,P=0,D=0,Q=0,S=-1,tol=sqrt(.Machine$double.eps)){ 
  n=length(data)
  constant=1:n   
  xmean=matrix(1,n,1) 
  if (d==0 & D==0) {
    fitit=arima(data, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S),
            xreg=xmean,include.mean=F, optim.control=list(trace=1,REPORT=1,reltol=tol))
} else if (xor(d==1, D==1)) {
    fitit=arima(data, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S),
            xreg=constant,include.mean=F, optim.control=list(trace=1,REPORT=1,reltol=tol))
} else fitit=arima(data, order=c(p,d,q), seasonal=list(order=c(P,D,Q), period=S), 
            optim.control=list(trace=1,REPORT=1,reltol=tol))
#
  if (S < 0) goof=20 else goof=3*S
  tsdiag(fitit,gof.lag=goof)
  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)
}