[Earlywarnings-commits] r5 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 14 12:59:23 CEST 2012
Author: vdakos
Date: 2012-04-14 12:59:23 +0200 (Sat, 14 Apr 2012)
New Revision: 5
Modified:
pkg/R/.Rhistory
pkg/R/generic_ews.R
Log:
update files
Modified: pkg/R/.Rhistory
===================================================================
--- pkg/R/.Rhistory 2012-01-17 16:18:03 UTC (rev 4)
+++ pkg/R/.Rhistory 2012-04-14 10:59:23 UTC (rev 5)
@@ -1,200 +1,947 @@
-defaults write org.R-project.R force.LANG en_US.UTF-8
-help("defaults")
-? corSpatial
-corSpatial
-? corSpatial
-logit
-x<-seq(1,12, by=0.1)
-xt<-logit(385.7*(x-5.41),min=0.91,max=0,94)
-library(Rcmdr)
-?logit
-x<-seq(1,12, by=0.1)
-xt<-logit(385.7*(x-5.41),min=0.91,max=0,94)
-xt<-logit(385.7*(x-5.41),min=0.91,max=0.94)
-xt<-logit(x,min=0.91,max=0.94)
-plot(xt)
-st
-xt
-xt<-logit(385.7*(x-5.41),min=0.91,max=0.94)
-xt
-xt<-logit(x,min=0,max=1)
-?logit
-x<-seq(0,10, by=0.1)
-xt<-logit(x,min=0,max=10)
-plot(xt)
-xt<-logit(385.7*(x-5.41),min=0.91,max=0.94)
-x<-seq(0.91,0.94, by=0.1)
-x<-seq(0.91,0.94, by=0.001)
-xt<-logit(385.7*(x-5.41),min=0.91,max=0.94)
-xt
-x<-seq(1,12, by=0.001)
-xt<-logit(385.7*(x-5.41),min=1,max=12)
-xt
-x<-seq(1,12, by=0.1)
-xt<-logit(385.7*(x-5.41),min=1,max=12)
-xt
-?logit
-xt<-logit(385.7*(x-5.41),min=0.91,max=0.94)
-xt
-source("/Users/vasilisdakos/Dropbox/current_projects/Paper_PhilTransSpIssue/R_trends_bandwidth_slidingwindow.r")
-rm(list=ls(all=TRUE))
+plot(xx)
+plot(xx[,1])
+Y=xx
+plot(Y)
+Y[1]
+Y[2]
+Y[1,]
+Y[1,1]
+plot(1:length(Y),Y)
+plot(ts(1:length(Y)),Y)
+length(Y)
+xx=read.table('CSD_6Dec11.txt',header=FALSE)
+length(xx)
+dim(xx)
+?length
+?data
+xx=data('CSD_6Dec11.txt')
+# Conditional Heteroskedasticity#
+# Author: Timothy Cline, October 25, 2011.#
+# Modified by: Vasilis Dakos, January 3, 2012.#
#
+ch_ews<-function(timeseries, winsize = 10, alpha=0.1, optim=TRUE,lags=4){#
+ #ts.in<-ts(timeseries) #strict data-types the input data as tseries object for use in later steps#
+ if (is.null(dim(timeseries)[2])){#
+ ts.in=timeseries#
+ timeindex=1:length(ts.in)#
+ }else if(dim(timeseries)[2]==2){#
+ ts.in=timeseries[,2]#
+ timeindex=timeseries[,1]#
+ }else{#
+ warning("not right format of timeseries input")#
+ }#
+ #
+ winSize=round(winsize*length(ts.in)/100)#
+ sto<-matrix(nrow=(length(ts.in)-(winSize-1)),ncol=5) # creates a matrix to store output#
+ #
+ count<-1 #place holder for writing to the matrix#
+ for(i2 in winSize:length(ts.in)){ # loop to iterate through the model values by window lengths of the input value#
+ #
+ #the next line applys the autoregressive model optimized using AIC #
+ #then we omit the first data point(s) which come back as NA #
+ #and square the residuals#
+ if(optim==TRUE){#
+ arm<-ar(ts.in[(i2-(winSize-1)):i2],method='ols')#
+ }else{#
+ arm<-ar(ts.in[(i2-(winSize-1)):i2],aic=FALSE,order.max=lags,method='ols') }#
+ resid1<-na.omit(arm$resid)^2#
+ #
+ l1<-length(resid1) # stores the number of residuals for many uses later#
+ lm1<-lm(resid1[2:l1]~resid1[1:(l1-1)])#calculates simple OLS model of describing t+1 by t#
+ #
+ # calculates the critical value: Chi-squared critical value using desired #
+ # alpha level and 1 degree of freedom / number of residuals used in regression#
+ critical<-qchisq((1-alpha),df=1)/(length(resid1)-1) #
+ #
+ sto[count,1]<-timeindex[i2] # stores a time component#
+ sto[count,2]<-summary(lm1)$r.squared # stores the r.squared for this window#
+ sto[count,3]<-critical # stores the critical value#
+ #
+ # the next flow control group stores a simple 1 for significant test or 0 for non-significant test#
+ if(summary(lm1)$r.squared>critical){#
+ sto[count,4]<-1#
+ }else{sto[count,4]<-0}#
+ sto[count,5]<-arm$order#
+ count<-count+1 # increment the place holder#
+ }#
+ #
+ sto<-data.frame(sto) # data types the matrix as a data frame#
+ colnames(sto)<-c("time","r.squared","critical.value","test.result","ar.fit.order") # applies column names to the data frame#
+ #
+ #This next series of flow control statements will adjust the max and minimum values to yield prettier plots#
+ #In some circumstances it is possible for all values to be far greater or far less than the critical value; in all cases we want the critical line ploted on the figure#
+ if(max(sto$r.squared)<critical){#
+ maxY<-critical+0.02 #
+ minY<-critical-0.02#
+ }else if(min(sto$r.squared>=critical)){#
+ minY<-critical-0.02#
+ maxY<-critical+0.02#
+ }else{#
+ maxY<-max(sto$r.squared)#
+ minY<-min(sto$r.squared)#
+ } #
+ #
+ #this creates a very simple plot that is well fitted to the data. #
+ #it also plots the critical value line #
+#
+ par(mar=(c(0,4,0,1)+0),oma=c(5,1,2,1),mfrow=c(2,1))#
+ plot(timeindex,ts.in,type="l",ylab="data",xlab="",cex.axis=0.8,cex.lab=0.8,xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ plot(timeindex[winSize:length(timeindex)],sto$r.squared,ylab=expression(paste("R^2")),xlab="time",type="b",cex=0.5,cex.lab=0.8, cex.axis=0.8,las=1,ylim=c(min(sto$r.squared),max(sto$r.squared)),xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","conditional heteroskedasticity",bty = "n")#
+ abline(h=sto$critical,lwd=0.5,lty=2,col=2)#
+ mtext("time",side=1,line=2,cex=0.8)#outer=TRUE print on the outer margin#
+ #
+ return(sto)#
+}
+head(xx)
+head(x)
+head(xxx)
+plot(xxx)
+ch_ews(xxx)
+ch_ews(ts(xxx))
+w=ts(xxx)
+plot(w)
+ch_ews(ts(xxx))
+ch_ews(w)
+dim(w)
+p=as.vector(w)
+head(p)
+ch_ews(p)
+?write
+write(p,'foldbif.txt')
+xx=read.delim('foldbif.txt',header=FALSE)
+dim(xx)
+head(xx)
+getwd()
+dim(p)
+length(p)
+write(p,'foldbif.txt',ncolumns=1)
+xx=read.delim('foldbif.txt',header=FALSE)
+head(xx)
+dim(xx)
+length(xx)
+plot(xx)
+plot(p)
+?read.table
+xx=read.delim('foldbif.txt',header=FALSE)
+plot(xx)
+xx=read.delim('foldbif.txt',header=TRUE)
+plot(xx)
+head(xx)
+xx=read.delim('foldbif.txt',header=FALSE)
+xx=scan('foldbif.txt',header=FALSE)
+xx=scan('foldbif.txt')
+plot(xx)
+dim(x)
+dim(xx)
+source("/Users/vasilisdakos/Dropbox/early_warning_tools_&_data_sets/R_code/Rpackage/earlywarnings/man/generic_ews.Rd")
+?plot
+?bds.test
+load("/Users/vasilisdakos/Documents/warningsignals/data/glaciationIII.rda")
+a=load("/Users/vasilisdakos/Documents/warningsignals/data/glaciationIII.rda")
+a
+a[1]
+a[2]
+a[2,]
+a[,2]
+?data
+?install.package
+?install.packages
+?data
+library(earlywarnings)
+?earlywarnings
+?generic_ews
+data("foldbif")
+foldbif
+plot(foldbif)
+plot(as.matrix(foldbif))
+data(circulation)
+plot(circulation)
+head(circulation)
+head(circulation)?
+?bandi5
+?ch_ews
+data(circulation)
+?generic_ews
+data(circulation)
+plot(circulation)
+x=data(circulation)
+dim(x)
+head()
+head(x)
+x
+head(circulation)
+generic_ews(circulation)
+lenght(circulation)
+length(circulation)
+dim(circlation)
+dim(circulation)
+x=circulation[1:782,]
+generic_ews(circulation)
+generic_ews(x)
+generic_ews(x,detrending="gaussian",ARn=FALSE)
+generic_ews(x,detrending="gaussian",ARn=FALSE)
+getwd()
+generic_ews(x,detrending="gaussian",ARn=FALSE)
+x
+getwd()
+dir
+list
+list()
+dir()
+generic_ews(x,detrending="gaussian",ARn=FALSE)
+# Generic Early Warning Signals#
+# Author: Vasilis Dakos, January 2, 2012#
+
+#
+generic_ews<-function(timeseries,winsize=50,detrending=c("no","gaussian","linear","first-diff"),bandwidth=NULL,logtransform=FALSE,interpolate=FALSE,ARn=FALSE,powerspectrum=FALSE){#
+ #
+ # Load required packages#
+ require('lmtest')#
+ require('nortest')#
+ require('stats')#
+ require('som')#
+ require('Kendall')#
+ require('KernSmooth')#
+ require('e1071')#
+ #
+ timeseries<-ts(timeseries) #strict data-types the input data as tseries object for use in later steps#
+ if (dim(timeseries)[2]==1){#
+ Y=timeseries#
+ timeindex=1:dim(timeseries)[1]#
+ }else if(dim(timeseries)[2]==2){#
+ Y<-timeseries[,2]#
+ timeindex<-timeseries[,1]#
+ }else{#
+ warning("not right format of timeseries input")#
+ }#
+ #
+ # Interpolation#
+ if (interpolate){#
+ YY<-approx(timeindex,Y,n=length(Y),method="linear")#
+ Y<-YY$y#
+ }else{#
+ Y<-Y}#
+ #
+ # Log-transformation#
+ if (logtransform){#
+ Y<-log(Y+1)}#
+ #
+ # Detrending #
+ detrending<-match.arg(detrending) #
+ if (detrending=="gaussian"){#
+ if (is.null(bandwidth)){#
+ bw<-round(bw.nrd0(timeindex))}else{#
+ bw<-round(length(Y)*bandwidth)/100}#
+ smYY<-ksmooth(timeindex,Y,kernel=c("normal"), bandwidth=bw, range.x=range(timeindex),n.points=length(timeindex))#
+ nsmY<-Y-smYY$y#
+ smY<-smYY$y#
+ }else if(detrending=="linear"){#
+ nsmY<-resid(lm(Y~timeindex))#
+ smY<-fitted(lm(Y~timeindex))#
+ }else if(detrending=="first-diff"){#
+ nsmY<-diff(Y)#
+ timeindexdiff<-timeindex[1:(length(timeindex)-1)]#
+ }else if(detrending=="no"){#
+ smY<-Y#
+ nsmY<-Y#
+ }
+#
#
-library(lmtest)
+
#
-library(nortest)
+ # Rearrange data for indicator calculation#
+ mw<-round(length(Y)*winsize)/100
#
-library(stats)
+ omw<-length(nsmY)-mw+1 ##number of moving windows
#
-library(som)
+ low<-6
#
-library(Kendall)
+ high<-omw
#
-library(KernSmooth)
+ nMR<-matrix(data=NA,nrow=mw,ncol=omw)
#
-library(e1071)
-library(lmtest)
-? acf
-?cor
-Y<-rand(500,1)
-?bds.test
-# Drift Diffusion Jump Nonparametrics Early Warning Signals#
-# Author: Stephen R Carpenter, 15 December 2011#
-# Modified by: Vasilis Dakos, January 4, 2012
+ x1<-1:mw
#
+ for (i in 1:omw){
+#
+ Ytw<-nsmY[i:(i+mw-1)]
+#
+ nMR[,i]<-Ytw}
+#
#
-require('KernSmooth')
+ # Calculate indicators
#
+ nARR<-numeric()
+#
+ nSD<-numeric()
+#
+ nSK<-numeric()
+#
+ nKURT<-numeric()
+#
+ nACF<-numeric()
+#
+ nDENSITYRATIO<-numeric()
+#
+ nSPECT<-matrix(0,nrow=omw, ncol=ncol(nMR))#
+ nCV<-numeric()#
+ smARall<-numeric()#
+ smARmaxeig<-numeric()#
+ detB<-numeric()#
+ ARn<-numeric()
+#
#
-# Function to compute Bandi, Johannes etc estimators for time series x
+ for (i in 1:ncol(nMR)){
#
-# Inputs:
+ nYR<-ar.ols(nMR[,i],aic= FALSE, order.max=1, dmean=FALSE, intercept=FALSE)
#
-# x0 is the regressor
+ nARR[i]<-nYR$ar
#
-# dx is the first difference of x0
+ nSD[i]<-sd(nMR[,i], na.rm = TRUE)
#
-# nx is number of first differences
+ nSK[i]<-abs(skewness(nMR[,i],na.rm=TRUE))
#
-# DT is time step
+ nKURT[i]<-kurtosis(nMR[,i],na.rm=TRUE)#
+ nCV[i]<-nSD[i]/mean(nMR[,i])
#
-# bw is the bandwidth for the kernel
+ ACF<-acf(nMR[,i], lag.max = 1, type = c("correlation"), plot=FALSE)
#
-# na is number of a values for computing the kernel
+ nACF[i]<-ACF$acf[2]
#
-# avec is the mesh for the kernel#
+ spectfft<-spec.ar(nMR[,i],n.freq=omw,plot=FALSE,order=1)
+#
+ nSPECT[,i]<-spectfft$spec
+#
+ nDENSITYRATIO[i]<- spectfft$spec[low]/spectfft$spec[high]#
+ #
+ ## RESILIENCE IVES 2003 Indicators based on AR(n)#
+ ARall<-ar.ols(nMR[,i],aic= TRUE,order.max=6,demean=F, intercept=F)#
+ smARall[i]<-ARall$ar[1]#
+ ARn[i]<-ARall$order#
+ roots<-Mod(polyroot(c(rev(-ARall$ar),1)))#
+ smARmaxeig[i]<-max(roots)#
+ detB[i]<-(prod(roots))^(2/ARn[i])}#
+ #
+ nRETURNRATE=1/nARR#
+
+#
+ # Estimate Kendall trend statistic for indicators
+#
+ timevec<-seq(1,length(nARR))
+#
+ KtAR<-cor.test(timevec,nARR,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtACF<-cor.test(timevec,nACF,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtSD<-cor.test(timevec,nSD,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtSK<-cor.test(timevec,nSK,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtKU<-cor.test(timevec,nKURT,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtDENSITYRATIO<-cor.test(timevec,nDENSITYRATIO,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)#
+ KtRETURNRATE<-cor.test(timevec,nRETURNRATE,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)#
+ KtCV<-cor.test(timevec,nCV,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
#
-Bandi5 <- function(x0,dx,nx,DT,bw,na,avec) {
+ # Plotting#
+ # Generic Early-Warnings
#
-# Set up constants and useful preliminaries
+ dev.new()#
+ par(mar=(c(0,2,0,1)+0),oma=c(7,2,3,1),mfrow=c(5,2))#
+ plot(timeindex,Y,type="l",ylab="",xlab="",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ if((detrending=="no") | (detrending=="first-diff")){#
+ }else{#
+ lines(timeindex,smY,col=2,xaxt="n",)}#
+ if(detrending=="no"){#
+ plot(timeindex,Y,ylab="",xlab="",xaxt="n",type="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ text(mean(timeindex),(max(Y)-min(Y))/2,"no residuals - no detrending")#
+ }else if (detrending=="first-diff"){#
+ limit<-max(c(max(abs(nsmY))))#
+ plot(timeindexdiff,nsmY,ylab="",xlab="",type="l",xaxt="n",las=1,ylim=c(- limit,limit),xlim=c(timeindexdiff[1],timeindexdiff[length(timeindexdiff)]))#
+ legend("topleft","first-differenced",bty = "n") }else{#
+ limit<-max(c(max(abs(nsmY))))#
+ plot(timeindex,nsmY,ylab="",xlab="",type="h",xaxt="n",las=1,ylim=c(- limit,limit),xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","residuals",bty = "n")}#
+ plot(timeindex[mw:length(nsmY)],nARR,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)])) #3#
+ legend("bottomleft",paste("Kendall tau=",round(KtAR$estimate,digits=3)),bty = "n")#
+ legend("topleft","ar(1)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nACF,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)])) #4#
+ legend("bottomleft",paste("Kendall tau=",round(KtACF$estimate,digits=3)),bty = "n")#
+ legend("topleft","acf(1)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nRETURNRATE,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtRETURNRATE$estimate,digits=3)),bty = "n")#
+ legend("topleft","return rate",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nDENSITYRATIO,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c (timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtDENSITYRATIO$estimate,digits=3)),bty = "n")#
+ legend("topleft","density ratio",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nSD,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtSD$estimate,digits=3)),bty = "n")#
+ legend("topleft","standard deviation",bty = "n")#
+ if(detrending=="no"){#
+ plot(timeindex[mw:length(nsmY)],nCV,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtCV$estimate,digits=3)),bty = "n")#
+ legend("topleft","coefficient of variation",bty = "n")}else{#
+ plot(timeindex,Y,ylab="",xlab="",type="n",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ text(mean(timeindex),mean(Y),"no coeff var estimated - data detrended")}#
+ plot(timeindex[mw:length(nsmY)],nSK,type="l",ylab="",xlab="",las=1,cex.lab=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","skewness",bty = "n")#
+ legend("bottomleft",paste("Kendall tau=",round(KtSK$estimate,digits=3)),bty = "n")#
+ mtext("time",side=1,line=2,cex=0.8)#
+ plot(timeindex[mw:length(nsmY)],nKURT,type="l",ylab="",xlab="",las=1,cex.lab=1,xlim=c(timeindex [1],timeindex[length(timeindex)]))#
+ legend("topleft","kurtosis",bty = "n")#
+ legend("bottomleft",paste("Kendall tau=",round(KtKU$estimate,digits=3)),bty = "n")#
+ mtext("time",side=1,line=2,cex=0.8)#
+ mtext("Generic Early-Warnings",side=3,line=0.2, outer=TRUE)#outer=TRUE print on the outer margin#
#
-SF <- 1/(bw*sqrt(2*pi)) # scale factor for kernel calculation
+ # Resilience Estimators based on AR(n)#
+ if (ARn==TRUE){#
+ dev.new()#
+ par(mar=(c(1,2,0,1)+0.2),oma=c(4,2,3,1),mfrow=c(2,2))#
+ plot(timeindex[mw:length(nsmY)],ARn,type="p",ylab="",xlab="",xaxt="n",cex=0.1,las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)])) #10#
+ legend("topleft","AR(n)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],smARmaxeig,type="l",ylab="",xlab="",xaxt="n",las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","max eigenvalue",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],detB,type="l",ylab="",xlab="",cex.lab=1,las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ mtext("time",side=1,line=2,cex=0.8)#
+ legend("topleft","geometric mean root AR(n)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],smARall,type="l",ylab="",xlab="",cex.lab=1,las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ mtext("time",side=1,line=2,cex=0.8)#
+ legend("topleft","b1 of AR(n)",bty = "n")
#
-x02 <- x0*x0 # second power of x
+ mtext("Resilience Estimators based on AR(n)",side=3,line=0.2, outer=TRUE)#
+ }#
+ #
+ # Power spectrum#
+ if (powerspectrum==TRUE){
#
-dx2 <- dx*dx # second power of dx
+ dev.new()#
+ par(mar=(c(4.6,4,0.5,2)+0.2),oma=c(0.5,1,2,1))
#
-dx4 <- dx2*dx2 # fourth power of dx
+ image(x=(spectfft$freq[2:length(spectfft$freq)]),y=(seq(1,ncol (nSPECT),by=1)),log(nSPECT[2:length(spectfft$freq),]),ylab="rolling window",xlab="frequency",log="x",xlim=c(spectfft$freq[2],spectfft$freq [length(spectfft$freq)]),col=topo.colors(20),xaxs="i")#
+ contour(x=(spectfft$freq[2:length(spectfft$freq)]),y=(seq(1,ncol (nSPECT),by=1)),log(nSPECT[2:length(spectfft$freq),]),add=TRUE)#
+ mtext("Power spectrum within rolling windows",side=3,line=0.2, outer=TRUE)#
+ }#
+ #
+ # Output#
+ out<-data.frame(timeindex[mw:length(nsmY)],nARR,nSD,nSK,nKURT,nCV,nRETURNRATE,nDENSITYRATIO,nACF)#
+ colnames(out)<-c("timeindex","ar1","sd","sk","kurt","cv","returnrate","densratio","acf1")#
+ return(out)#
+ }
+generic_ews(x,detrending="gaussian",ARn=FALSE)
+generic_ews(x)
+# Generic Early Warning Signals#
+# Author: Vasilis Dakos, January 2, 2012#
+
#
-dx6 <- dx2*dx4 # sixth power of dx
+generic_ews<-function(timeseries,winsize=50,detrending=c("no","gaussian","linear","first-diff"),bandwidth=NULL,logtransform=FALSE,interpolate=FALSE,ARn=FALSE,powerspectrum=FALSE){#
+ #
+ # Load required packages#
+ require('lmtest')#
+ require('nortest')#
+ require('stats')#
+ require('som')#
+ require('Kendall')#
+ require('KernSmooth')#
+ require('e1071')#
+ #
+ timeseries<-ts(timeseries) #strict data-types the input data as tseries object for use in later steps#
+ if (dim(timeseries)[2]==1){#
+ Y=timeseries#
+ timeindex=1:dim(timeseries)[1]#
+ }else if(dim(timeseries)[2]==2){#
+ Y<-timeseries[,2]#
+ timeindex<-timeseries[,1]#
+ }else{#
+ warning("not right format of timeseries input")#
+ }#
+ #
+ # Interpolation#
+ if (interpolate){#
+ YY<-approx(timeindex,Y,n=length(Y),method="linear")#
+ Y<-YY$y#
+ }else{#
+ Y<-Y}#
+ #
+ # Log-transformation#
+ if (logtransform){#
+ Y<-log(Y+1)}#
+ #
+ # Detrending #
+ detrending<-match.arg(detrending) #
+ if (detrending=="gaussian"){#
+ if (is.null(bandwidth)){#
+ bw<-round(bw.nrd0(timeindex))}else{#
+ bw<-round(length(Y)*bandwidth)/100}#
+ smYY<-ksmooth(timeindex,Y,kernel=c("normal"), bandwidth=bw, range.x=range(timeindex),n.points=length(timeindex))#
+ nsmY<-Y-smYY$y#
+ smY<-smYY$y#
+ }else if(detrending=="linear"){#
+ nsmY<-resid(lm(Y~timeindex))#
+ smY<-fitted(lm(Y~timeindex))#
+ }else if(detrending=="first-diff"){#
+ nsmY<-diff(Y)#
+ timeindexdiff<-timeindex[1:(length(timeindex)-1)]#
+ }else if(detrending=="no"){#
+ smY<-Y#
+ nsmY<-Y#
+ }
#
-# Compute matrix of kernel values
+
#
-Kmat <- matrix(0,nrow=na,ncol=nx)
+
#
-for(i in 1:(nx)) { # loop over columns (x0 values)
+ # Rearrange data for indicator calculation#
+ mw<-round(length(Y)*winsize)/100
#
- Kmat[,i] <- SF*exp(-0.5*(x0[i]-avec)*(x0[i]-avec)/(bw*bw))
+ omw<-length(nsmY)-mw+1 ##number of moving windows
#
- }
+ low<-6
#
-# Compute M1, M2, M4, moment ratio and components of variance for each value of a
+ high<-omw
#
-M1.a <- rep(0,na)
+ nMR<-matrix(data=NA,nrow=mw,ncol=omw)
#
-M2.a <- rep(0,na)
+ x1<-1:mw
#
-M4.a <- rep(0,na)
+ for (i in 1:omw){
#
-M6M4r <- rep(0,na) # vector to hold column kernel-weighted moment ratio
+ Ytw<-nsmY[i:(i+mw-1)]
#
-mean.a <- rep(0,na) # centering of conditional variance
+ nMR[,i]<-Ytw}
#
-SS.a <- rep(0,na) # sum of squares
+
#
-for(i in 1:na) { # loop over rows (a values)
+ # Calculate indicators
#
- Ksum <- sum(Kmat[i,]) # sum of weights
+ nARR<-numeric()
#
- M1.a[i] <- (1/DT)*sum(Kmat[i,]*dx)/Ksum
+ nSD<-numeric()
#
- M2.a[i] <- (1/DT)*sum(Kmat[i,]*dx2)/Ksum
+ nSK<-numeric()
#
- M4.a[i] <- (1/DT)*sum(Kmat[i,]*dx4)/Ksum
+ nKURT<-numeric()
#
- M6.c <- (1/DT)*sum(Kmat[i,]*dx6)/Ksum
+ nACF<-numeric()
#
- M6M4r[i] <- M6.c/M4.a[i]
+ nDENSITYRATIO<-numeric()
#
- mean.a[i] <- sum(Kmat[i,]*x0[2:(nx+1)])/Ksum
+ nSPECT<-matrix(0,nrow=omw, ncol=ncol(nMR))#
+ nCV<-numeric()#
+ smARall<-numeric()#
+ smARmaxeig<-numeric()#
+ detB<-numeric()#
+ ARn<-numeric()
#
- SS.a[i] <- sum(Kmat[i,]*x02[2:(nx+1)])/Ksum
+
#
- }
+ for (i in 1:ncol(nMR)){
#
-# Compute conditional variance
+ nYR<-ar.ols(nMR[,i],aic= FALSE, order.max=1, dmean=FALSE, intercept=FALSE)
#
-S2.x <- SS.a - (mean.a*mean.a) # sum of squares minus squared mean
+ nARR[i]<-nYR$ar
#
-# Compute jump frequency, diffusion and drift
+ nSD[i]<-sd(nMR[,i], na.rm = TRUE)
#
-sigma2.Z <- mean(M6M4r)/(5) # average the column moment ratios
+ nSK[i]<-abs(skewness(nMR[,i],na.rm=TRUE))
#
-lamda.Z <- M4.a/(3*sigma2.Z*sigma2.Z)
+ nKURT[i]<-kurtosis(nMR[,i],na.rm=TRUE)#
+ nCV[i]<-nSD[i]/mean(nMR[,i])
#
-sigma2.dx <- M2.a - (lamda.Z*sigma2.Z)
+ ACF<-acf(nMR[,i], lag.max = 1, type = c("correlation"), plot=FALSE)
#
-# set negative diffusion estimates to zero
+ nACF[i]<-ACF$acf[2]
#
-diff.a <- ifelse(sigma2.dx>0,sigma2.dx,0)
+ spectfft<-spec.ar(nMR[,i],n.freq=omw,plot=FALSE,order=1)
#
-sigma2.dx <- M2.a # total variance of dx
+ nSPECT[,i]<-spectfft$spec
#
-mu.a <- M1.a
+ nDENSITYRATIO[i]<- spectfft$spec[low]/spectfft$spec[high]#
+ #
+ ## RESILIENCE IVES 2003 Indicators based on AR(n)#
+ ARall<-ar.ols(nMR[,i],aic= TRUE,order.max=6,demean=F, intercept=F)#
+ smARall[i]<-ARall$ar[1]#
+ ARn[i]<-ARall$order#
+ roots<-Mod(polyroot(c(rev(-ARall$ar),1)))#
+ smARmaxeig[i]<-max(roots)#
+ detB[i]<-(prod(roots))^(2/ARn[i])}#
+ #
+ nRETURNRATE=1/nARR#
+
#
-outlist <- list(mu.a,sigma2.dx,diff.a,sigma2.Z,lamda.Z,S2.x)
+ # Estimate Kendall trend statistic for indicators
#
-# outputs of function:
+ timevec<-seq(1,length(nARR))
#
-# mu.a is drift
+ KtAR<-cor.test(timevec,nARR,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
#
-# sigma2.dx is total variance of dx
+ KtACF<-cor.test(timevec,nACF,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
#
-# diff.a is diffusion
+ KtSD<-cor.test(timevec,nSD,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
#
-# sigma2.Z is jump magnitude
+ KtSK<-cor.test(timevec,nSK,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
#
-# lamda.Z is jump frequency
+ KtKU<-cor.test(timevec,nKURT,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
#
-# S2.x is conditional variance
+ KtDENSITYRATIO<-cor.test(timevec,nDENSITYRATIO,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)#
+ KtRETURNRATE<-cor.test(timevec,nRETURNRATE,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)#
+ KtCV<-cor.test(timevec,nCV,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
#
-return(outlist)
+
#
-} # end Bandi function
+ # Plotting#
+ # Generic Early-Warnings
#
+ dev.new()#
+ par(mar=(c(0,2,0,1)+0),oma=c(7,2,3,1),mfrow=c(5,2))#
+ plot(timeindex,Y,type="l",ylab="",xlab="",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ if((detrending=="no") | (detrending=="first-diff")){#
+ }else{#
+ lines(timeindex,smY,col=2,xaxt="n",)}#
+ if(detrending=="no"){#
+ plot(timeindex,Y,ylab="",xlab="",xaxt="n",type="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ text(mean(timeindex),(max(Y)-min(Y))/2,"no residuals - no detrending")#
+ }else if (detrending=="first-diff"){#
+ limit<-max(c(max(abs(nsmY))))#
+ plot(timeindexdiff,nsmY,ylab="",xlab="",type="l",xaxt="n",las=1,ylim=c(- limit,limit),xlim=c(timeindexdiff[1],timeindexdiff[length(timeindexdiff)]))#
+ legend("topleft","first-differenced",bty = "n") }else{#
+ limit<-max(c(max(abs(nsmY))))#
+ plot(timeindex,nsmY,ylab="",xlab="",type="h",xaxt="n",las=1,ylim=c(- limit,limit),xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","residuals",bty = "n")}#
+ plot(timeindex[mw:length(nsmY)],nARR,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)])) #3#
+ legend("bottomleft",paste("Kendall tau=",round(KtAR$estimate,digits=3)),bty = "n")#
+ legend("topleft","ar(1)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nACF,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)])) #4#
+ legend("bottomleft",paste("Kendall tau=",round(KtACF$estimate,digits=3)),bty = "n")#
+ legend("topleft","acf(1)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nRETURNRATE,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtRETURNRATE$estimate,digits=3)),bty = "n")#
+ legend("topleft","return rate",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nDENSITYRATIO,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c (timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtDENSITYRATIO$estimate,digits=3)),bty = "n")#
+ legend("topleft","density ratio",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],nSD,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtSD$estimate,digits=3)),bty = "n")#
+ legend("topleft","standard deviation",bty = "n")#
+ if(detrending=="no"){#
+ plot(timeindex[mw:length(nsmY)],nCV,ylab="",xlab="",type="l",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("bottomleft",paste("Kendall tau=",round(KtCV$estimate,digits=3)),bty = "n")#
+ legend("topleft","coefficient of variation",bty = "n")}else{#
+ plot(timeindex,Y,ylab="",xlab="",type="n",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ text(mean(timeindex),mean(Y),"no coeff var estimated - data detrended")}#
+ plot(timeindex[mw:length(nsmY)],nSK,type="l",ylab="",xlab="",las=1,cex.lab=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","skewness",bty = "n")#
+ legend("bottomleft",paste("Kendall tau=",round(KtSK$estimate,digits=3)),bty = "n")#
+ mtext("time",side=1,line=2,cex=0.8)#
+ plot(timeindex[mw:length(nsmY)],nKURT,type="l",ylab="",xlab="",las=1,cex.lab=1,xlim=c(timeindex [1],timeindex[length(timeindex)]))#
+ legend("topleft","kurtosis",bty = "n")#
+ legend("bottomleft",paste("Kendall tau=",round(KtKU$estimate,digits=3)),bty = "n")#
+ mtext("time",side=1,line=2,cex=0.8)#
+ mtext("Generic Early-Warnings",side=3,line=0.2, outer=TRUE)#outer=TRUE print on the outer margin#
+#
+ # Resilience Estimators based on AR(n)#
+ if (ARn){#
+ dev.new()#
+ par(mar=(c(1,2,0,1)+0.2),oma=c(4,2,3,1),mfrow=c(2,2))#
+ plot(timeindex[mw:length(nsmY)],ARn,type="p",ylab="",xlab="",xaxt="n",cex=0.1,las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)])) #10#
+ legend("topleft","AR(n)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],smARmaxeig,type="l",ylab="",xlab="",xaxt="n",las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ legend("topleft","max eigenvalue",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],detB,type="l",ylab="",xlab="",cex.lab=1,las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ mtext("time",side=1,line=2,cex=0.8)#
+ legend("topleft","geometric mean root AR(n)",bty = "n")#
+ plot(timeindex[mw:length(nsmY)],smARall,type="l",ylab="",xlab="",cex.lab=1,las=1,cex.axis=0.8,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ mtext("time",side=1,line=2,cex=0.8)#
+ legend("topleft","b1 of AR(n)",bty = "n")
+#
+ mtext("Resilience Estimators based on AR(n)",side=3,line=0.2, outer=TRUE)#
+ }#
+ #
+ # Power spectrum#
+ if (powerspectrum){
+#
+ dev.new()#
+ par(mar=(c(4.6,4,0.5,2)+0.2),oma=c(0.5,1,2,1))
+#
+ image(x=(spectfft$freq[2:length(spectfft$freq)]),y=(seq(1,ncol (nSPECT),by=1)),log(nSPECT[2:length(spectfft$freq),]),ylab="rolling window",xlab="frequency",log="x",xlim=c(spectfft$freq[2],spectfft$freq [length(spectfft$freq)]),col=topo.colors(20),xaxs="i")#
+ contour(x=(spectfft$freq[2:length(spectfft$freq)]),y=(seq(1,ncol (nSPECT),by=1)),log(nSPECT[2:length(spectfft$freq),]),add=TRUE)#
+ mtext("Power spectrum within rolling windows",side=3,line=0.2, outer=TRUE)#
+ }#
+ #
+ # Output#
+ out<-data.frame(timeindex[mw:length(nsmY)],nARR,nSD,nSK,nKURT,nCV,nRETURNRATE,nDENSITYRATIO,nACF)#
+ colnames(out)<-c("timeindex","ar1","sd","sk","kurt","cv","returnrate","densratio","acf1")#
+ return(out)#
+ }
+generic_ews(x)
+# Generic Early Warning Signals#
+# Author: Vasilis Dakos, January 2, 2012#
+
+#
+generic_ews<-function(timeseries,winsize=50,detrending=c("no","gaussian","linear","first-diff"),bandwidth=NULL,logtransform=FALSE,interpolate=FALSE,AR_n=FALSE,powerspectrum=FALSE){#
+ #
+ # Load required packages#
+ require('lmtest')#
+ require('nortest')#
+ require('stats')#
+ require('som')#
+ require('Kendall')#
+ require('KernSmooth')#
+ require('e1071')#
+ #
+ timeseries<-ts(timeseries) #strict data-types the input data as tseries object for use in later steps#
+ if (dim(timeseries)[2]==1){#
+ Y=timeseries#
+ timeindex=1:dim(timeseries)[1]#
+ }else if(dim(timeseries)[2]==2){#
+ Y<-timeseries[,2]#
+ timeindex<-timeseries[,1]#
+ }else{#
+ warning("not right format of timeseries input")#
+ }#
+ #
+ # Interpolation#
+ if (interpolate){#
+ YY<-approx(timeindex,Y,n=length(Y),method="linear")#
+ Y<-YY$y#
+ }else{#
+ Y<-Y}#
+ #
+ # Log-transformation#
+ if (logtransform){#
+ Y<-log(Y+1)}#
+ #
+ # Detrending #
+ detrending<-match.arg(detrending) #
+ if (detrending=="gaussian"){#
+ if (is.null(bandwidth)){#
+ bw<-round(bw.nrd0(timeindex))}else{#
+ bw<-round(length(Y)*bandwidth)/100}#
+ smYY<-ksmooth(timeindex,Y,kernel=c("normal"), bandwidth=bw, range.x=range(timeindex),n.points=length(timeindex))#
+ nsmY<-Y-smYY$y#
+ smY<-smYY$y#
+ }else if(detrending=="linear"){#
+ nsmY<-resid(lm(Y~timeindex))#
+ smY<-fitted(lm(Y~timeindex))#
+ }else if(detrending=="first-diff"){#
+ nsmY<-diff(Y)#
+ timeindexdiff<-timeindex[1:(length(timeindex)-1)]#
+ }else if(detrending=="no"){#
+ smY<-Y#
+ nsmY<-Y#
+ }
+#
#
-# MAIN FUNCTION#
-ddjnonparam_ews<-function(timeseries,bandwidth=0.6,na=500,logtransform=TRUE,interpolate=FALSE){#
+
#
- # Create timeindex vector and data vector#
- if (is.null(dim(timeseries)[2])){#
- Y<-timeseries#
- timeindex<-1:length(timeseries)#
+ # Rearrange data for indicator calculation#
+ mw<-round(length(Y)*winsize)/100
+#
+ omw<-length(nsmY)-mw+1 ##number of moving windows
+#
+ low<-6
+#
+ high<-omw
+#
+ nMR<-matrix(data=NA,nrow=mw,ncol=omw)
+#
+ x1<-1:mw
+#
+ for (i in 1:omw){
+#
+ Ytw<-nsmY[i:(i+mw-1)]
+#
+ nMR[,i]<-Ytw}
+#
+
+#
+ # Calculate indicators
+#
+ nARR<-numeric()
+#
+ nSD<-numeric()
+#
+ nSK<-numeric()
+#
+ nKURT<-numeric()
+#
+ nACF<-numeric()
+#
+ nDENSITYRATIO<-numeric()
+#
+ nSPECT<-matrix(0,nrow=omw, ncol=ncol(nMR))#
+ nCV<-numeric()#
+ smARall<-numeric()#
+ smARmaxeig<-numeric()#
+ detB<-numeric()#
+ ARn<-numeric()
+#
+
+#
+ for (i in 1:ncol(nMR)){
+#
+ nYR<-ar.ols(nMR[,i],aic= FALSE, order.max=1, dmean=FALSE, intercept=FALSE)
+#
+ nARR[i]<-nYR$ar
+#
+ nSD[i]<-sd(nMR[,i], na.rm = TRUE)
+#
+ nSK[i]<-abs(skewness(nMR[,i],na.rm=TRUE))
+#
+ nKURT[i]<-kurtosis(nMR[,i],na.rm=TRUE)#
+ nCV[i]<-nSD[i]/mean(nMR[,i])
+#
+ ACF<-acf(nMR[,i], lag.max = 1, type = c("correlation"), plot=FALSE)
+#
+ nACF[i]<-ACF$acf[2]
+#
+ spectfft<-spec.ar(nMR[,i],n.freq=omw,plot=FALSE,order=1)
+#
+ nSPECT[,i]<-spectfft$spec
+#
+ nDENSITYRATIO[i]<- spectfft$spec[low]/spectfft$spec[high]#
+ #
+ ## RESILIENCE IVES 2003 Indicators based on AR(n)#
+ ARall<-ar.ols(nMR[,i],aic= TRUE,order.max=6,demean=F, intercept=F)#
+ smARall[i]<-ARall$ar[1]#
+ ARn[i]<-ARall$order#
+ roots<-Mod(polyroot(c(rev(-ARall$ar),1)))#
+ smARmaxeig[i]<-max(roots)#
+ detB[i]<-(prod(roots))^(2/ARn[i])}#
+ #
+ nRETURNRATE=1/nARR#
+
+#
+ # Estimate Kendall trend statistic for indicators
+#
+ timevec<-seq(1,length(nARR))
+#
+ KtAR<-cor.test(timevec,nARR,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtACF<-cor.test(timevec,nACF,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtSD<-cor.test(timevec,nSD,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtSK<-cor.test(timevec,nSK,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtKU<-cor.test(timevec,nKURT,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+ KtDENSITYRATIO<-cor.test(timevec,nDENSITYRATIO,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)#
+ KtRETURNRATE<-cor.test(timevec,nRETURNRATE,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)#
+ KtCV<-cor.test(timevec,nCV,alternative=c("two.sided"),method=c("kendall"),conf.level=0.95)
+#
+
+#
+ # Plotting#
+ # Generic Early-Warnings
+#
+ dev.new()#
+ par(mar=(c(0,2,0,1)+0),oma=c(7,2,3,1),mfrow=c(5,2))#
+ plot(timeindex,Y,type="l",ylab="",xlab="",xaxt="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ if((detrending=="no") | (detrending=="first-diff")){#
+ }else{#
+ lines(timeindex,smY,col=2,xaxt="n",)}#
+ if(detrending=="no"){#
+ plot(timeindex,Y,ylab="",xlab="",xaxt="n",type="n",las=1,xlim=c(timeindex[1],timeindex[length(timeindex)]))#
+ text(mean(timeindex),(max(Y)-min(Y))/2,"no residuals - no detrending")#
+ }else if (detrending=="first-diff"){#
+ limit<-max(c(max(abs(nsmY))))#
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/earlywarnings -r 5
More information about the Earlywarnings-commits
mailing list