tperiod <- function(x, figure=FALSE)

#        COMPUTING PERIODOGRAM AND DETECTING SIGNIFICANT PERIODICITIES
#                      USING FISHER'S AND SIEGEL'S TEST
#        Additive model: x(t)=y(t)+e(t) is assumed where
#        y(t) ... deterministic component with periodicities to be detected
#        e(t) ... error component, e: IID(0,sigma2)
# =============================================================================
# INPUT:
# x   ... is an 1 x n vector of observations
# fig ... logical
#         figure handle for the periodogram with emphasized significant
#         harmonics
#         FALSE=default=no plot
# OUTPUT:
# list containing the following components:
# I ... periodogram: vector of length n2
#                    n2=trunc((n-1)/2)
# kF ... indices of significant harmonics using Fisher's test
# kS ... indices of significant harmonics using Siegel's test

# Iveta Hellebrandova, 2011
# according to the MATLAB function "tperiod":
#   (C) 1997-1999, V. Vesely, Masaryk University of Brno (Czech Republic)
#   Rev. 1, Dec. 22, 1999

{
if (min(dim(as.matrix(x)))>1) stop('x is not a vector')
n<-length(x)
m<-trunc((n-1)/2)              # length of the periodogram
if (m<5) stop('x must be longer than 10')
m<-min(m,1000)
I<-fft(x)
I<-(2/n)*abs(I[2:(m+1)])^2    # Computing the periodogram
II<-I
lF<-numeric(); lS<-numeric(); k<-numeric()
# source('fisher.r')
while (m>4) {                  # loop to detect significant periodicities
   f <- fisher(m)              # Fisher's and Siegel's 95# percentile
   gf <- f$pfisher
   tsi <- f$psiegel
   Y <- II/sum(II)
   W <- max(Y)                # Fisher's statistic
   kmax <- which.max(Y)
   Ts<-sum(pmax(0, Y - 0.6*gf))    # Siegel's statistic
   if (length(lF)==0 & W<=gf) lF <- length(k)
   if (length(lS)==0 & Ts<=tsi) lS <- length(k)
   if (W>gf | Ts>tsi) {
      k <- c(k,kmax[1])
      II[kmax[1]] <- 0     # scratch component just detected
      m <- m-1
    } else {
    break
    }
  }

if (figure) {
  plot(I)
  lines(I, type='h', lty='dashed')
  points(k,I[k],pch=8)
  text(k,I[k],round(n/k,1), pos=4)
  abline(h=0)
  if (lS>lF)
   title("PERIODOGRAM (Siegel's test)")  else
   title("PERIODOGRAM (Fisher's test)")
  }
kF <- k[1:lF]
kS <- k[1:lS]
return(list(I=I, kF=kF, kS=kS))
}
#========================================================================
fisher <- function(m)
#    95# PERCENTILES FOR THE FISHER'S AND SIEGEL'S TEST OF PERIODICITY
# =======================================================================
# INPUT:
# m ... vector of lengths of the periodogram, 5<=m[j]<=1000
# OUTPUT:
# list containing the following components:
# pfisher... pfisher[j]=Fisher's 95% percentile for m[j]
# psiegel... psiegel[j]=Siegel's 95% percentile for m[j]

# Iveta Hellebrandova, 2011
# according to the MATLAB function "fisher":
#   (C) 1997-1999, V. Vesely, Masaryk University of Brno (Czech Republic)
#   Rev. 1, Dec. 22, 1999

{
if (any(m<5)) stop('Periodogram must have length at least 5')
if (any(m>1000)) stop('The periodogram length must not exceed 1000')
GF <- rbind(c(0.6837722340, 0.274),
      c(0.4449525606, 0.1812916218),
      c(0.3346119552, 0.1401751436),
      c(0.2704042202, 0.1159434232),
      c(0.2280456085, 0.0997325432),
      c(0.1978406246, 0.0880198043),
      c(0.1751294860, 0.0791060851),
      c(0.1573827449, 0.0720634801),
      c(0.1431035999, 0.0663391976),
      c(0.1313472721, 0.0615819252),
      c(0.1130887896, 0.0541010364),
      c(0.0995257984, 0.0484581309),
      c(0.0890258722, 0.0440297152),
      c(0.0806402063, 0.0404493552),
      c(0.0737781851, 0.0374866267),
      c(0.0522136336, 0.0279346469),
      c(0.0407367370, 0.0226517317),
      c(0.0335542444, 0.0192460210),
      c(0.0286124993, 0.0168450268),
      c(0.0249931178, 0.0150496264),
      c(0.0222217969, 0.0136497353),
      c(0.0200280827, 0.0125235715),
      c(0.0182461407, 0.0115953645),
      c(0.0155221997, 0.0101493846),
      c(0.0135334433, 0.0090692380),
      c(0.0120144724, 0.0082276189),
      c(0.0108145470, 0.0075508836),
      c(0.0098415283,  0.0069932836))

mGF <- c(seq(5,50,5),seq(60,100,10),seq(150,500,50),seq(600,1000,100))
pfisher <- spline(mGF,GF[,1],xout=m)$y
psiegel <- spline(mGF,GF[,2],xout=m)$y
return(list(pfisher=pfisher, psiegel=psiegel))
}
#=======================================================================
plot.band <- function (x,lwd=1.5,type="o",pch=20, ...) {
  plot(x, type=type,pch=pch,lwd=lwd, ...)
  a <- time(x)
  i1 <- floor(min(a))
  i2 <- ceiling(max(a))
  y1 <- par('usr')[3]
  y2 <- par('usr')[4]
  if( par("ylog") ){
    y1 <- 10^y1
    y2 <- 10^y2
  }
  for (i in seq(from=i1, to=i2-1, by=2)) {
    polygon( c(i,i+1,i+1,i),
             c(y1,y1,y2,y2),
             col = 'grey',
             border = NA )
  }
  par(new=T)
  plot(x, lwd=lwd,type=type,pch=pch,...)
}
#=========================================================================
damsleth<-function(X,pstlevel=0.05)
# iterativni metoda nalezeni skrytych period pomoci
# nelinearni regrese (Gauss-Newton metoda)
# pocatecni hodnoty neznamych parametru se urcuji pomoci
# Fisherova testu periodogramu

{
if (min(dim(as.matrix(X)))>1) stop('X is not a vector')
n<-length(X);
TIME<-1:n
m<-as.integer(n/2)
Y<-X
aa<-NULL;bb<-NULL;freq<-NULL
konec<-FALSE
krok<-1
while (!konec & krok<=m)
{
   # krok 1
   # pomoci Fisherova testu se urci prvotni priblizeni
   beta0<-fishper(Y,pstlevel)

   # krok 2
   # nelinearni regrese
   if (is.null(beta0))
      konec<-TRUE
   else {
      modelnls<-nls(Y~trigpol1(TIME,A,B,freq),
                      start=list(A=beta0$A,B=beta0$B,freq=beta0$freq))
      aa<-c(aa,coef(modelnls)[1])
      bb<-c(bb,coef(modelnls)[2])
      freq<-c(freq,coef(modelnls)[3])
      Y<-resid(modelnls)
      krok<-krok+1
   }
}
return(list(A=aa,B=bb,freq=freq))
}
#=======================================================================
fishper<-function(XX,pstlevel=0.05)
# funkce pro urceni nejvice signifikantni periody-Fisheruv test
# vypocet periodogramu
{
# vytvoreni omega(k)=(2*pi/n)*k  k-1,...,n/2

if (min(dim(as.matrix(XX)))>1) stop('X is not a vector')
n<-length(XX)
if (n %% 2==0)
   {X<-XX[-1]
   n<-n-1
   } else
   X<-XX

m<-as.integer(n/2)

omega<-matrix((1:m),ncol=1)*2*pi/n
#                 n
# a(k)=sqrt(2/n)*suma X(t)*cos(omega(k)*t))
#                 t=1
a<-sqrt(2/n)*cos(kronecker(matrix((1:n),nrow=1),omega))%*%X
#                 n
# b(k)=sqrt(2/n)*suma X(t)*sin(omega(k)*t))
#                 t=1
b<-sqrt(2/n)*sin(kronecker(matrix((1:n),nrow=1),omega))%*%X
I<-1/(4*pi)*(a^2+b^2)
a<-a*sqrt(2/n)
b<-b*sqrt(2/n)
# Fisheruv test pro maximalni hodnotu periodogramu
Vo<-max(I)
ind<-which.max(I)
sumV<-sum(I)
pstG<-0
Wo<-Vo/sumV
# Pst(W>Wo)=m*(1-1*Wo)^(m-1)+(m nad 2)*(1-2*Wo)^(m-1)+....
# secita se tak dlouho,dokud (1-k*Wo)>0
k<-1
pstG<-0;
sign<-1
while (k*Wo<1)
   {pstG<-pstG+sign*comnum(m,k)*(1-k*Wo)^(m-1)
   k<-k+1
   sign<-sign*(-1)
   }
ifelse(pstG<=pstlevel,
         return(list(A=a[ind],B=b[ind],freq=omega[ind])),
         return(NULL))
}
#============================================================
comnum<-function(n,k)
# kombinacni cislo (n nad k) = (n/1)*((n-1)/2)*...*((n-k+1)/k)
comb<-prod(seq(n,n-k+1,by=-1)/(1:k))
#=============================================================
trigpol1<-function(TIME,A,B,freq)
# funkce pro vypocet trigonometrickeho polynomu s jednou frekvenci
# beta0=list(A, B, freq)
# y(t)=A*cos(TIME*freq)+B*sin(TIME*freq)
{
y<-A*cos(freq*TIME)+B*sin(freq*TIME)
}
#=============================================================
trigoval<-function(TIME,A,B,freq)
#function Y=trigoval(T,aa,bb,freq)
# trigonometricky polynom z vyznamnych period
#           p
# trig(t)= sum(aa(i)*cos((freq(i)*t)+bb(i)*sin((freq(i)*t)
#          i=1
{
if (min(dim(as.matrix(TIME)))>1) stop('TIME is not a vector')
Y<-cos(kronecker(matrix(TIME,ncol=1),
                 matrix(freq,nrow=1)))%*%matrix(A,ncol=1)+
   sin(kronecker(matrix(TIME,ncol=1),
                 matrix(freq,nrow=1)))%*%matrix(B,ncol=1)
}
#=================================================================
trigofit<-function(XX,pstlevel=0.05)
#function [I,omega,aa,bb,freq,amplituda,faze,L,a,b]=trigofit(XX,pstlevel)
# funkce pro nalezeni  amplitud a fazi skrytych period
# vypocet periodogramu
# vytvoreni omega(k)=(2*pi/n)*k  k-1,...,n/2
{
if (min(dim(as.matrix(XX)))>1) stop('X is not a vector')
n<-length(XX)
if (n %% 2==0)
   {X<-XX[-1]
   n<-n-1
   } else
   X<-XX
m<-as.integer(n/2)
omega<-matrix((1:m),ncol=1)*2*pi/n
L<-rep(FALSE,m)
a<-sqrt(2/n)*cos(kronecker(matrix((1:n),nrow=1),omega))%*%X
b<-sqrt(2/n)*sin(kronecker(matrix((1:n),nrow=1),omega))%*%X
I<-1/(4*pi)*(a^2+b^2)
a<-a*sqrt(2/n)
b<-b*sqrt(2/n)
amplituda=sqrt(a^2+b^2)*sqrt(2/n)
faze=tan(b/a)

# Metoda skrytych period
o<-order(I)
V<-I[o]
sumV<-sum(V)
i<-m
pstG<-0
pocet<-0
aa<-NULL
bb<-NULL
freq<-NULL
while ((pstG<=pstlevel)&(i<=m))
   {
   Wo<-V[i]/sumV
   sumV<-sumV-V[i]
   # Pst(W>Wo)=m*(1-1*Wo)^(m-1)+(m nad 2)*(1-2*Wo)^(m-1)+....
   # secita se tak dlouho,dokud (1-k*Wo)>0
   k<-1
   pstG<-0
   sign<-1
   while (k*Wo<1)
   {
      pstG<-pstG+sign*comnum(m,k)*(1-k*Wo)^(m-1);
      k<-k+1;
      sign<-sign*(-1);
   }
   if (pstG<=pstlevel)
      {
      L[o[i]]=TRUE
      aa<-c(aa,a[o[i]])
      bb<-c(bb,b[o[i]])
      freq=c(freq,omega[o[i]])
      pocet<-pocet+1
      }
   i<-i-1
   }
return(list(periodogram=as.vector(I),freq=as.vector(omega),
            signA=aa,signB=bb,signFreq=freq,
            amplituda=as.vector(amplituda),faze=as.vector(faze),signL=L,
            A=as.vector(a),B=as.vector(b)))
}
#=======================================================================
fndptr<-
function(y,x=(1:length(y))-0.5*(1+length(y)),alpha=0.05,maxdgr=10,mindgr=1)
{
x<-x-mean(x)
vh<-NULL;vsk<-NULL;vak<-NULL
vaick<-NULL;vsrk<-NULL;vhqk2<-NULL;vhqk3<-NULL
ny<-length(y)
maxdgr<-min(maxdgr,ny-2)
mindgr<-max(mindgr,1)
for(k in mindgr:maxdgr)
{
 if(k==1) pomf<-as.formula("y~x") else
          pomf<-as.formula(paste("y~x",paste("+I(x^",2:k,")",sep="",collapse=""),sep=""))
 m<-lm(pomf)
 CI<-confint(m)[k+1,]
 if((0>=CI[1])&(0<=CI[2])) hyp<-0 else hyp<-1
 sk1<-sum(m$residuals^2)/(ny-k-1)
   vh<-c(vh,hyp)
   vsk<-c(vsk,sk1)
   vak<-c(vak,sk1*(1+k*(ny^(-0.25))))
   vaick<-c(vaick,log(sk1)+(2*k/ny))
   vsrk<-c(vsrk,log(sk1)+k*log(ny)/ny)
   vhqk2<-c(vhqk2,log(sk1)+4*k*log(log(ny))/ny)
   vhqk3<-c(vhqk3,log(sk1)+6*k*log(log(ny))/ny)
}
vd<-mindgr:maxdgr
L=vh>0
vysl<-data.frame(vd,vh,vsk,vak,vaick,vsrk,vhqk2,vhqk3)
names(vysl)<-c("dgr","significant","S_k","A_k","AIC_k","SR_k","HQ_k(c=2)","HQ_k(c=3)")

TXT<-c("S_k (Mean Square Error)","A_k","AIC_k","SR_k","HQ_k(c=2)","HQ_k(c=3)")

par(mfrow=c(2,3),mar = c(2, 2, 1, 0) + 0.5)
for(kk in 1:6)
{
plot(vysl[,1],vysl[,kk+2],type="o",main=TXT[kk],xlab="k",ylab="")
if(kk==1 & sum(L)>0) text(vysl[L,1],vysl[L,kk+2],
               labels=vysl[L,1],pos=3)
ind<-which.min(vysl[,kk+2])
abline(v=vysl[ind,1],lty=2,col="red")
points(vysl[ind,1],vysl[ind,kk+2],pch=19,cex=1.25)
mtext(paste("opt =",vysl[ind,1]),cex=0.85,line=-1)
}
return(vysl)
}
#===================================================================================
powtr<-
function(x,seglen=8,figure=FALSE,location="median",variability="iqr",
         figure2=FALSE,figure3=FALSE)
#       COMPUTING ESTIMATE OF PARAMETER lambda FOR THE POWER TRANSFORM:
#                        y = x.^lambda   if lambda~=0
#                          = log(x)      if lambda =0  (natural logarithm)
# ===================================================================================
# INPUT:
# x      ... is an n x 1 or 1 x n data vector of POSITIVE observed data
#            or vector of smoothed values (determnistic component) provided
#            that the residual error component r describing the decomposition
#            x = y + r is supplied instead of seglen
# seglen ... segment length used for computation of "std. deviation" growth
#            (recommended: 4<=seglen<=12)
#             default segment length=8
# fig ... figure handle for the plot of estimated "std. deviation" versus mean
#         optional: []=default=no plot
# OUTPUT:
# lambdaE ... =[lambda,E] is and 1 x 2 vector, where
#              lambda is estimate of parameter lambda
#              [lambda-E,lambda+E] is its confidence interval of 95#
# y  ... =data transformed according to parameter lambda
#         vector of the same size as x

# (C) 2011, M. Forbelska, Masaryk University of Brno (Czech Republic)
#

{
if (min(dim(as.matrix(x)))>1) stop('x is not a vector')
n<-length(x)

if(seglen < 4)
   {
      warning("segment length too short: 4 will be used")
      seglen<-4
   }
if(seglen > 12)
   {
      warning("segment length too long: 12 will be used")
      seglen<-12
   }
#---------------------------
   nseg <- floor(n/seglen)
   if(nseg < 5)
     {
      warning("less than 5 segments: minimal segment length 4 will be used")
      seglen <- 4
      nseg <- floor(n./seglen)
      if(nseg < 5)
         {
          stop("Still less than 5 segments - estimate impossible")
         }
     }
   ns <- nseg*seglen
   n0 <- n - ns
   gr<-rep(1:nseg,each=seglen)
   if(n0>4) app<-rep(nseg+1,n0) else app<-rep(nseg,n0)
   gr<-c(gr,app)
   MINx<-min(x)
   if(MINx<=0)
      {
         delta<-0.005*diff(range(x))-MINx
         x<-x+delta
         MSG<-paste("Minimum data value <= 0 so  -min(x)+0.005*diff(range(x))=",delta,
                    "  added to all values",sep="")

         message(MSG)
       }

   if(location == "mean") LOC<-tapply(x,gr,mean) else LOC<-tapply(x,gr,median)

   if(variability == "sd") VAR<-tapply(x,gr,sd) else
      if(variability == "range") VAR<-tapply(x,gr,function(x) diff(range(x)))  else
         VAR<-tapply(x,gr,IQR)

   # pro normalni rozdeleni IQR = 1.34896 sigma

   data<-data.frame(logLocation=log(LOC),logVariability=log(VAR))
   lm.fit<-lm(logVariability~logLocation,data=data)
   #A<-coef(lm.fit)[1]
   B<-summary(lm.fit)$coefficients[2,]
   lambda<-1-B[1]
   lambdaL<-lambda-qnorm(0.975)*B[2]
   lambdaH<-lambda+qnorm(0.975)*B[2]

   if(1>lambdaL & 1lambdaL & 0>>",TXT,"<<<"))

   if (figure)
      {
       rLoc<-with(data,range(logLocation))
       nn<-100
       new <- data.frame(logLocation = seq(rLoc[1],rLoc[2],length.out=nn))
       pred.w.plim <- predict(lm.fit, new, interval="prediction")
       pred.w.clim <- predict(lm.fit, new, interval="confidence")
       par(mar=c(4,4,1,0)+0.5)
       with(data,plot(logVariability~logLocation))
       matlines(new$logLocation,cbind(pred.w.clim, pred.w.plim[,-1]),col=c(2,3,3,4,4),
                lty=c(1,2,2,3,3), type="l")
       mtext(paste("b = ",round(B[1],6),",  lambda = ",round(lambda,6),
                   ",  CI = (",round(lambdaL,6)," , ",round(lambdaH,6),")",sep=""),side=1,line=-1)
       mtext(paste("seglen = ",seglen,",  location = ",location,",   variability = ",variability,sep=""),line=-1)
       title(main=TXT,cex.main=01)
     }

   if (figure2)
      {
      xTS<-ts(x,frequency=seglen)
      #dev.new()
       par(mar=c(2,4,1,0)+0.5,mfrow=c(2,1))
      plot.band(xTS,lwd=1,type="p",pch=1,col="gray30",ylab="original values")
      for(ik in 1:(nseg))
       {
        polygon(c(ik,ik+1,ik+1,ik),c(rep(LOC[ik]-0.5*VAR[ik],2),rep(LOC[ik]+0.5*VAR[ik],2)),col="yellow",border="yellow")
        lines(c(ik,ik+1),rep(LOC[ik],2),col="red",lwd=2.5)
       }
      if(n0>4)
        {
        polygon(c(ik+1,ik+1+n0/seglen,ik+1+n0/seglen,ik+1),
                c(rep(LOC[ik+1]-0.5*VAR[ik+1],2),rep(LOC[ik+1]+0.5*VAR[ik+1],2)),
                col="yellow",border="yellow")
        lines(c(ik+1,ik+1+n0/seglen),rep(LOC[ik+1],2),col="red",lwd=2.5)
        } else { if(n0>0)
                    {
                     polygon(c(ik+1,ik+1+n0/seglen,ik+1+n0/seglen,ik+1),
                             c(rep(LOC[ik]-0.5*VAR[ik],2),rep(LOC[ik]+0.5*VAR[ik],2)),
                             col="yellow",border="yellow")
                     lines(c(ik+1,ik+1+n0/seglen),rep(LOC[ik],2),col="red",lwd=2.5)
                     }
                }
      points(time(xTS),xTS,pch=1)
       mtext(paste("seglen = ",seglen,",  location = ",location,",   variability = ",variability,sep=""))
     #---------
      x<-transfx
      xTS<-ts(x,frequency=seglen)
      if(location == "mean") LOC<-tapply(x,gr,mean) else LOC<-tapply(x,gr,median)
      if(variability == "sd") VAR<-tapply(x,gr,sd) else
        if(variability == "range") VAR<-tapply(x,gr,function(x) diff(range(x)))  else  VAR<-tapply(x,gr,IQR)
      plot.band(xTS,lwd=1,type="p",pch=1,col="gray30",ylab="transformed values")
      for(ik in 1:(nseg))
       {
        polygon(c(ik,ik+1,ik+1,ik),c(rep(LOC[ik]-0.5*VAR[ik],2),rep(LOC[ik]+0.5*VAR[ik],2)),col="yellow",border="yellow")
        lines(c(ik,ik+1),rep(LOC[ik],2),col="red",lwd=2.5)
       }
      if(n0>4)
        {
        polygon(c(ik+1,ik+1+n0/seglen,ik+1+n0/seglen,ik+1),
                c(rep(LOC[ik+1]-0.5*VAR[ik+1],2),rep(LOC[ik+1]+0.5*VAR[ik+1],2)),
                col="yellow",border="yellow")
        lines(c(ik+1,ik+1+n0/seglen),rep(LOC[ik+1],2),col="red",lwd=2.5)
        } else { if(n0>0)
                    {
                     polygon(c(ik+1,ik+1+n0/seglen,ik+1+n0/seglen,ik+1),
                             c(rep(LOC[ik]-0.5*VAR[ik],2),rep(LOC[ik]+0.5*VAR[ik],2)),
                             col="yellow",border="yellow")
                     lines(c(ik+1,ik+1+n0/seglen),rep(LOC[ik],2),col="red",lwd=2.5)
                     }
                }
      points(time(xTS),xTS,pch=1)
       mtext(paste("lambda = ",round(lambda,6),
                   ",  CI = (",round(lambdaL,6)," , ",round(lambdaH,6),")",sep=""),side=1,line=-1)
       mtext(TXT)
      }
#---------
   if (figure3)
      {
       par(mar=c(2,4,1,0)+0.5,mfrow=c(2,1))
       plot(x~factor(gr))
              mtext(paste("seglen = ",seglen,",  location = ",location,",   variability = ",variability,sep=""))

       plot(transfx~factor(gr))
       mtext(paste("lambda = ",round(lambda,6),
                   ",  CI = (",round(lambdaL,6)," , ",round(lambdaH,6),")",sep=""),side=1,line=-1)
       mtext(TXT)
      }
#---------
return(list(lambda=lambda,transfx=transfx,txt=TXT))
}
#==============================================================================================================
boxplotSegments<-
function(x,seglen=8,shift=0,appendL=0,appendH=0,...)
# seglen    delka segmentu (alespon 4)
# shift     posunuti pri deleni x na segmenty
# appendL ... ma smysl pro shift>0
#   2   nechat samostane
#   0   vynechat
#   1   spojit s prvni skupinou delky seglen
# appendH ... zbytek < seglen
#   2   nechat samostane
#   0   vynechat
#   1   spojit s posledni skupinou delky seglen
{
if (min(dim(as.matrix(x)))>1) stop('x is not a vector')
x<-as.vector(x)
n<-length(x)
if(seglen < 4)
   {
      warning("segment length too short: 4 will be used")
      seglen<-4
   }
#------------------------------
nseg <- floor((n-shift)/seglen)
ns <- nseg*seglen
n0 <- n - ns - shift
#------------------------------
idStart<-shift+1
idEnd<-n-n0
grL<-NULL
grH<-NULL
igr<-1
#------------------------------
if(shift>0 & appendL==2)
   {
    idStart<-1
    nseg<-nseg+1
    igr<-2
    grL<-rep(1,shift)
   }
if(shift>0 & appendL==1)
   {
    idStart<-1
    grL<-rep(1,shift)
    }
if(n0>0 & appendH==2)
   {
    idEnd<-n
    grH<-rep(nseg+1,n0)
   }
if(n0>0 & appendH==1)
   {
    idEnd<-n
    grH<-rep(nseg,n0)
   }
gr<-c(grL,rep(igr:nseg,each=seglen),grH)
X<-x[idStart:idEnd]
segments<-factor(gr)
#-------------------------
plot(X~segments,...)
}
#============================================================================================================
HistFit<-
function(x,main="Histogram, Kernel Density Estimate, Normal Curve",
         colH="bisque",ticksize=0.02,colKD="red",colN="dodgerblue",
         ltyKD=2,ltyN=1,lwdKD=2,lwdN=2,positLegend=NULL,...)
{
if (min(dim(as.matrix(x)))>1) stop('x is not a vector')
x<-as.vector(x)
#---
if(is.null(positLegend))
{
xx <- x[!is.na(x)]
nn <- length(xx)
Skewness<-(sum((xx - mean(xx))^3)/nn)/(sum((xx - mean(xx))^2)/nn)^(3/2)
if(Skewness < 0) positLegend<-"topleft" else positLegend<-"topright"
}
#---
xhist <- hist(x, breaks="Sturges", plot=FALSE)
nHx<-length(xhist$breaks)
xr<-c(xhist$breaks[1],xhist$breaks[nHx])
nD<-300
xdens<-density(x,n=nD,from=xr[1],to=xr[2])
nxfit<-dnorm(xdens$x,mean=mean(x),sd=sd(x))
topx <- max(c(xhist$density,nxfit,xdens$y))
#---
par(mar=c(4,2,1,0)+0.75)
h<-hist(x,probability=TRUE, breaks="Sturges",
        ylim=c(0,topx),main=main,col=colH,...)
lines(xdens$x,nxfit, col=colN,lty=ltyN,lwd=lwdN)
lines(xdens$x,xdens$y,lwd=lwdKD,col=colKD,lty=ltyKD)
rug(x,side=1,col="grey20",ticksize=ticksize)
#---
legend(positLegend,legend=c("Histogram","Kernel Density Estimate","Normal Density"),
       col=c("black",colKD,colN),lty=c(1,ltyKD,ltyN),bty="n")
}
#============================================================================================================
CenterFilter<-function(k)
  {
  if(k%%2!=0) stop("k must be a odd number")
  centerFilter<-c(1/(2*k),rep(1/k,k-1),1/(2*k))
  }
#============================================================================================================
 SzSmallTrendOld<-function(xts,figure=TRUE)
 {
 if(class(xts) != "ts") stop("xts must be a time series (ts)")
 d<-frequency(xts)
 n<-length(xts)
 m<-ceiling(n/d)
 shift<-start(xts)[2]-1
 r<-m*d-n-shift
 Season<-factor(as.integer(gl(d,1,n))-shift)
 Year<-factor(floor(time(xts)))
 data<-data.frame(x=as.vector(xts),grS=Season,grY=Year)
 y<-as.vector(xts)
 dd<-diag(rep(1,d))[,1:(d-1)]
 dd[d,]<-rep(-1,d-1)
 pm<-matrix(rep(1,m),ncol=1)
 Xd<-kronecker(pm,dd)[(1+shift):(n+shift),]
 pd<-matrix(rep(1,d),ncol=1)
 dm<-diag(rep(1,m))
 Xm<-kronecker(dm,pd)[(1+shift):(n+shift),]
 X<-cbind(Xm,Xd)
 print(summary(M1<-lm(y~X-1)))
 M1.mj<-dm%*%coef(M1)[1:m]
 M1.sk<-dd%*%coef(M1)[(m+1):(m+d-1)]
 if(figure)
 {
 tt<-as.vector(time(xts))
 xx<-as.vector(xts)
 fit<-fitted(M1)
 tr<-rep(M1.mj,each=d,length.out=n)
 ylim<-range(range(xx),range(fit),tr)
 par(mfrow=c(1,1),mar=c(2,2,2,0)+0.5)
 plot(tt,xx,type="o",pch=20,col="black",ylim=ylim)
 lines(tt,fit,type="o",pch=22,col="red",bg="yellow",cex=0.85)
 lines(tt,tr,type="s",col="dodgerblue",lwd=2)
 abline(v=as.numeric(as.character(levels(data$grY))),col="gray",lty=2)
 legend(par("usr")[1],par("usr")[4],bty="n",
       legend=c("original","estimate"),lty=c(1,1),
       pch=c(20,22),bg=c("black","yellow"),col=c("black","red"))
  }
return(list(model=M1,mj=M1.mj,sk=M1.sk,fit=fit,
            mjts=rep(M1.mj,each=d,length.out=n),skts=rep(M1.sk,m,length.out=n)))
 }
#===============================================================================
 SzSmallTrendModif<-function(xts,figure=TRUE)
 {
 if(class(xts) != "ts") stop("xts must be a time series (ts)")
 d<-frequency(xts)
 n<-length(xts)
 m<-ceiling(n/d)
 shift<-start(xts)[2]-1
 r<-m*d-n-shift
 Season<-factor(as.integer(gl(d,1,n))-shift)
 Year<-factor(floor(time(xts)))
 data<-data.frame(x=as.vector(xts),grS=Season,grY=Year)   
 M1m<-lm(x~grY+grS,data,contrasts=list(grY=contr.sum,grS=contr.sum))
 M1m.terms<-predict(M1m,type="terms")
 M1m.sk<-M1m.terms[1:d,2]
 M1m.mj<-attr(M1m.terms,"constant")+M1m.terms[seq(1,n,d),1]
 if(figure)
 {
 tt<-as.vector(time(xts))
 xx<-as.vector(xts)

 fit<-fitted(M1m)
  tr<-rep(M1m.mj,each=d,length.out=n)
 ylim<-range(range(xx),range(fit),tr)
 par(mfrow=c(1,1),mar=c(2,2,2,0)+0.5)
 plot(tt,xx,type="o",pch=20,col="black",ylim=ylim)
 lines(tt,fit,type="o",pch=22,col="red",bg="yellow",cex=0.85)
 lines(tt,tr,type="s",col="dodgerblue",lwd=2)
 abline(v=as.numeric(as.character(levels(data$grY))),col="gray",lty=2)
 abline(h=attr(M1m.terms,"constant"),lty=2,col="gray")
 legend(FindPositionForLegend(xx),bty="n",
       legend=c("original","estimate"),lty=c(1,1),
       pch=c(20,22),bg=c("black","yellow"),col=c("black","red"))
  }
return(list(model=M1m,mu=coef(M1m)[1],mj=M1m.mj,sk=M1m.sk))
 }
#===============================================================================
FindPositionForLegend<-function(x)
{
if (min(dim(as.matrix(x)))>1) stop('x is not a vector')
n<-length(x)
m<-trunc(n/3)
MaxT<-c(max(x[1:m]),max(x[(m+1):(n-m)]),max(x[(n-m+1):n]))
MinB<-c(min(x[1:m]),min(x[(m+1):(n-m)]),min(x[(n-m+1):n]))
indT<-which.min(MaxT)
indB<-which.max(MinB)
if(max(abs(diff(MaxT[c(1:3,1)]))) > max(abs(diff(MinB[c(1:3,1)]))))
   PositionLegend<-c("topleft","top","topright")[indT]  else
   PositionLegend<-c("bottomleft","bottom","bottomright")[indB]
return(PositionLegend)
}
#===============================================================================
 SzSmallTrend<-function(xts,figure=TRUE)
 {
 if(class(xts) != "ts") stop("xts must be a time series (ts)")
 d<-frequency(xts)
 n<-length(xts)
 m<-ceiling(n/d)
 shift<-start(xts)[2]-1
 r<-m*d-n-shift
 Season<-factor(as.integer(gl(d,1,n))-shift)
 Year<-factor(floor(time(xts)))
 data<-data.frame(x=as.vector(xts),grS=Season,grY=Year)
 M11<-lm(x~grY+grS-1,data,contrasts=list(grY=diag(1,length(levels(data$grY))),grS=contr.sum))
 M11.terms<-predict(M11,type="terms")
 M11.sk<-M11.terms[1:d,2]
 M11.mj<-M11.terms[seq(1,n,d),1]
 if(figure)
 {
 tt<-as.vector(time(xts))
 xx<-as.vector(xts)
 fit<-fitted(M11)
 tr<-rep(M11.mj,each=d,length.out=n)
 ylim<-range(range(xx),range(fit),tr)
 par(mfrow=c(1,1),mar=c(2,2,2,0)+0.5)
 plot(tt,xx,type="o",pch=20,col="black",ylim=ylim)
 lines(tt,fit,type="o",pch=22,col="red",bg="yellow",cex=0.85)
 lines(tt,tr,type="s",col="dodgerblue",lwd=2)
 abline(v=as.numeric(as.character(levels(data$grY))),col="gray",lty=2)
 legend(FindPositionForLegend(xx),bty="n",
       legend=c("original","estimate"),lty=c(1,1),
       pch=c(20,22),bg=c("black","yellow"),col=c("black","red"))
  }
return(list(model=M11,mj=M11.mj,sk=M11.sk))
 }
#===============================================================================
SzPolTrend<-function(xts,Dgr=1,figure=TRUE)
{
if(class(xts) != "ts") stop("xts must be a time series (ts)")
nn<-300
d<-frequency(xts)
n<-length(xts)
m<-ceiling(n/d)
shift<-start(xts)[2]-1
r<-m*d-n-shift
grS<-factor(as.integer(gl(d,1,n))-shift)
grY<-factor(floor(time(xts)))
Time<-(1:n)-mean(1:n)
Data<-data.frame(x=as.vector(xts),grS=grS,Time=Time)
T.terms<-"Time"
if(Dgr>1) T.terms<-c(T.terms,paste("I(Time^",2:Dgr,")",sep=""))
strF<-paste("x~",paste(T.terms,collapse="+"),"+grS",sep="")
pomf<-as.formula(strF)
M<-lm(pomf,Data,contrasts=list(grS=contr.sum))
m.terms.grS<-predict(M,type="terms",terms="grS")
m.sk<-m.terms.grS[1:d]
if(figure)
{
tt<-as.vector(time(xts))
xx<-as.vector(xts)
TTime<-Time
TT<-tt
if(np.min ) {
    par(mfrow=c(5,1))
  } else {
    par(mfrow=c(4,1))
  }
  if(!is.ts(x))
    x <- ts(x)
  plot(x, axes=FALSE);
  axis(2); axis(3); box(lwd=2)
  if(bands) {
    a <- time(x)
    i1 <- floor(min(a))
    i2 <- ceiling(max(a))
    y1 <- par('usr')[3]
    y2 <- par('usr')[4]
    if( par("ylog") ){
      y1 <- 10^y1
      y2 <- 10^y2
    }
    for (i in seq(from=i1, to=i2-1, by=2)) {
      polygon( c(i,i+1,i+1,i), c(y1,y1,y2,y2),
               col='grey', border=NA )
    }
    lines(x)
  }
  acf(x, axes=FALSE)
  axis(2, las=2)
  box(lwd=2)
  mtext("ACF", side=2, line=2.5)
  pacf(x, axes=FALSE)
  axis(2, las=2)
  box(lwd=2)
  mtext("PACF", side=2, line=2.5)
  spectrum(x, col=par('fg'), log="dB",
           main="", axes=FALSE )
  axis(2, las=2)
  box(lwd=2)
  mtext("Spectrum", side=2, line=2.5)
  abline(v=1, lty=2, lwd=2)
  abline(v=2:10, lty=3)
  abline(v=1/2:5, lty=3)
  if( max(p)>p.min ) {
    main <-
    plot(p, type='h', ylim=c(0,1),
         lwd=3, main="", axes=F)
    axis(2, las=2)
    box(lwd=2)
    mtext("Ljung-Box p-value", side=2, line=2.5)
    abline(h=c(0,.05),lty=3)
  }
  par(op)
}