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) 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="FD", 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="FD", 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(n