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)
}