[Earlywarnings-commits] r3 - pkg pkg/R pkg/data pkg/man www www/images
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jan 17 17:12:59 CET 2012
Author: vdakos
Date: 2012-01-17 17:12:59 +0100 (Tue, 17 Jan 2012)
New Revision: 3
Added:
pkg/DESCRIPTION
pkg/R/
pkg/R/.Rhistory
pkg/R/bdstest_ews.R
pkg/R/ch_ews.R
pkg/R/ddjnonparam_ews.R
pkg/R/earlywarnings-internal.R
pkg/R/generic_ews.R
pkg/R/sensitivity_ews.R
pkg/R/surrogates_ews.R
pkg/data/
pkg/data/.Rhistory
pkg/data/YD2PB_grayscale.txt
pkg/data/circulation.txt
pkg/data/foldbif.txt
pkg/man/
pkg/man/.Rhistory
pkg/man/bdstest_ews.Rd
pkg/man/ch_ews.Rd
pkg/man/ddjnonparam_ews.Rd
pkg/man/earlywarnings-package.Rd
pkg/man/generic_ews.Rd
pkg/man/sensitivity_ews.Rd
pkg/man/surrogates_ews.Rd
www/images/
www/images/fig_website_R.pdf
Log:
add files
Added: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION (rev 0)
+++ pkg/DESCRIPTION 2012-01-17 16:12:59 UTC (rev 3)
@@ -0,0 +1,10 @@
+Package: earlywarnings
+Type: Package
+Title: Early Warning Signals Toolbox for Detecting Critical Transitions in Timeseries
+Version: 1.0
+Date: 2012-01-02
+Author: Vasilis Dakos, with contributions from S.C. Carpenter, T. Cline, A. Ives, V. Livina
+Maintainer: Vasilis Dakos <vasilis.dakos at gmail.com>
+Description: The Early-Warning-Signals Toolbox provides methods for estimating statistical changes in timeseries that can be used for identifying nearby critical transitions. Based on Dakos et al (2012)
+URL: http://www.early-warning-signals.org, http://www.vasilisdakos.net
+License: GPL (>= 2)
\ No newline at end of file
Added: pkg/R/.Rhistory
===================================================================
--- pkg/R/.Rhistory (rev 0)
+++ pkg/R/.Rhistory 2012-01-17 16:12:59 UTC (rev 3)
@@ -0,0 +1,558 @@
+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))
+#
+
+#
+library(lmtest)
+#
+library(nortest)
+#
+library(stats)
+#
+library(som)
+#
+library(Kendall)
+#
+library(KernSmooth)
+#
+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
+#
+
+#
+require('KernSmooth')
+#
+
+#
+# Function to compute Bandi, Johannes etc estimators for time series x
+#
+# Inputs:
+#
+# x0 is the regressor
+#
+# dx is the first difference of x0
+#
+# nx is number of first differences
+#
+# DT is time step
+#
+# bw is the bandwidth for the kernel
+#
+# na is number of a values for computing the kernel
+#
+# avec is the mesh for the kernel#
+
+#
+Bandi5 <- function(x0,dx,nx,DT,bw,na,avec) {
+#
+# Set up constants and useful preliminaries
+#
+SF <- 1/(bw*sqrt(2*pi)) # scale factor for kernel calculation
+#
+x02 <- x0*x0 # second power of x
+#
+dx2 <- dx*dx # second power of dx
+#
+dx4 <- dx2*dx2 # fourth power of dx
+#
+dx6 <- dx2*dx4 # sixth power of dx
+#
+# Compute matrix of kernel values
+#
+Kmat <- matrix(0,nrow=na,ncol=nx)
+#
+for(i in 1:(nx)) { # loop over columns (x0 values)
+#
+ Kmat[,i] <- SF*exp(-0.5*(x0[i]-avec)*(x0[i]-avec)/(bw*bw))
+#
+ }
+#
+# Compute M1, M2, M4, moment ratio and components of variance for each value of a
+#
+M1.a <- rep(0,na)
+#
+M2.a <- rep(0,na)
+#
+M4.a <- rep(0,na)
+#
+M6M4r <- rep(0,na) # vector to hold column kernel-weighted moment ratio
+#
+mean.a <- rep(0,na) # centering of conditional variance
+#
+SS.a <- rep(0,na) # sum of squares
+#
+for(i in 1:na) { # loop over rows (a values)
+#
+ Ksum <- sum(Kmat[i,]) # sum of weights
+#
+ M1.a[i] <- (1/DT)*sum(Kmat[i,]*dx)/Ksum
+#
+ M2.a[i] <- (1/DT)*sum(Kmat[i,]*dx2)/Ksum
+#
+ M4.a[i] <- (1/DT)*sum(Kmat[i,]*dx4)/Ksum
+#
+ M6.c <- (1/DT)*sum(Kmat[i,]*dx6)/Ksum
+#
+ M6M4r[i] <- M6.c/M4.a[i]
+#
+ mean.a[i] <- sum(Kmat[i,]*x0[2:(nx+1)])/Ksum
+#
+ SS.a[i] <- sum(Kmat[i,]*x02[2:(nx+1)])/Ksum
+#
+ }
+#
+# Compute conditional variance
+#
+S2.x <- SS.a - (mean.a*mean.a) # sum of squares minus squared mean
+#
+# Compute jump frequency, diffusion and drift
+#
+sigma2.Z <- mean(M6M4r)/(5) # average the column moment ratios
+#
+lamda.Z <- M4.a/(3*sigma2.Z*sigma2.Z)
+#
+sigma2.dx <- M2.a - (lamda.Z*sigma2.Z)
+#
+# set negative diffusion estimates to zero
+#
+diff.a <- ifelse(sigma2.dx>0,sigma2.dx,0)
+#
+sigma2.dx <- M2.a # total variance of dx
+#
+mu.a <- M1.a
+#
+outlist <- list(mu.a,sigma2.dx,diff.a,sigma2.Z,lamda.Z,S2.x)
+#
+# outputs of function:
+#
+# mu.a is drift
+#
+# sigma2.dx is total variance of dx
+#
+# diff.a is diffusion
+#
+# sigma2.Z is jump magnitude
+#
+# lamda.Z is jump frequency
+#
+# S2.x is conditional variance
+#
+return(outlist)
+#
+} # end Bandi function
+#
+
+#
+# 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)#
+ }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)}#
+
+#
+ # Preliminaries#
+ Xvec1<-Y#
+ Tvec1<-timeindex
+#
+ dXvec1 <- diff(Y)
+#
+
+#
+ DT <- Tvec1[2]-Tvec1[1]
+#
+ bw <- bandwidth*sd(Xvec1) # bandwidth
+#
+ alow <- min(Xvec1)
+#
+ ahigh <- max(Xvec1)
+#
+ na <- na
+#
+ avec <- seq(alow,ahigh,length.out=na)
+#
+ nx <- length(dXvec1)
+#
+
+#
+ # Bandi-type estimates
+#
+ ParEst <- Bandi5(Xvec1,dXvec1,nx,DT,bw,na,avec)
+#
+ Drift.vec <- ParEst[[1]]
+#
+ TotVar.dx.vec <- ParEst[[2]]
+#
+ Diff2.vec <- ParEst[[3]]
+#
+ Sigma2Z <- ParEst[[4]]
+#
+ LamdaZ.vec <- ParEst[[5]]
+#
+ S2.vec <- ParEst[[6]]
+#
+
+#
+ # Interpolate time courses of indicators
+#
+ TotVar.i <- approx(x=avec,y=TotVar.dx.vec,xout=Xvec1)
+#
+ TotVar.t <- TotVar.i$y
+#
+ Diff2.i <- approx(x=avec,y=Diff2.vec,xout=Xvec1)
+#
+ Diff2.t <- Diff2.i$y
+#
+ Lamda.i <- approx(x=avec,y=LamdaZ.vec,xout=Xvec1)
+#
+ Lamda.t <- Lamda.i$y
+#
+ S2.i <- approx(x=avec,y=S2.vec,xout=Xvec1)
+#
+ S2.t <- S2.i$y
+#
+
+#
+ # Plot the data
+#
+ dev.new()
+#
+ par(mfrow=c(2,1),mar=c(3, 3, 2, 2),mgp=c(1.5,0.5,0),oma=c(1,1,1,1))
+#
+ plot(Tvec1,Xvec1,type='l',col='black',lwd=2,xlab='',ylab='original data')
+#
+ grid()#
+ plot(Tvec1[1:length(Tvec1)-1],dXvec1,type='l',col='black',lwd=2,xlab='time',ylab='first-diff data')
+#
+ grid()#
+
+#
+ # Plot indicators versus a
+#
+ dev.new()
+#
+ par(mfrow=c(2,2),mar=c(3, 3, 2, 2) ,cex.axis=1,cex.lab=1,mgp=c(2,1,0),oma=c(1,1,2,1))
+#
+ plot(avec,S2.vec,type='l',lwd=1,col='black',xlab='a',ylab='conditional variance')
+#
+ plot(avec,TotVar.dx.vec,type='l',lwd=1,col='blue',xlab='a',ylab='total variance of dx')
+#
+ plot(avec,Diff2.vec,type='l',lwd=1,col='green',xlab='a',ylab='diffusion')
+#
+ plot(avec,LamdaZ.vec,type='l',lwd=1,col='red',xlab='a',ylab='jump intensity')
+#
+ mtext("DDJ nonparametrics versus a",side=3,line=0.1,outer=TRUE)#
+
+#
+ # Plot indicators versus time
+#
+ dev.new()
+#
+ par(mfrow=c(2,2),mar=c(3, 3, 2, 2),cex.axis=1,cex.lab=1,mgp=c(1.5,0.5,0),oma=c(1,1,2,1))
+#
+ plot(Tvec1,S2.t,type='l',lwd=1,col='black',xlab='time',ylab='conditional variance')
+#
+ plot(Tvec1,TotVar.t,type='l',lwd=1,col='blue',xlab='time',ylab='total variance of dx')
+#
+ plot(Tvec1,Diff2.t,type='l',lwd=1,col='green',xlab='time',ylab='diffusion')
+#
+ plot(Tvec1,Lamda.t,type='l',lwd=1,col='red',xlab='time',ylab='jump intensity')#
+ mtext("DDJ nonparametrics versus time",side=3,line=0.1,outer=TRUE)#
+#
+ # Output#
+ nonpar_x<-data.frame(S2.vec,TotVar.dx.vec,Diff2.vec,LamdaZ.vec)#
+ #out$x<-nonpar_x#
+#
+ nonpar_t<-data.frame(S2.t,TotVar.t,Diff2.t,Lamda.t)#
+ #out$t<-nonpar_t#
+ return(data.frame(nonpar_x,nonpar_t))#
+ }
+?bds.test
+require(tseries)
+?bds.test
+Y=rnorm(100)
+Y
+Y=rnorm(500)
+ss=ddjnonparam_ews(timeseries,bandwidth=0.6,na=500,logtransform=TRUE,interpolate=FALSE)
+ss=ddjnonparam_ews(Y,bandwidth=0.6,na=500,logtransform=TRUE,interpolate=FALSE)
+Y=rnorm(500)+10
+ss=ddjnonparam_ews(Y,bandwidth=0.6,na=500,logtransform=TRUE,interpolate=FALSE)
+ss
+head(ss)
+?data
+getwd()
+xx=read.table('CSD_6Dec11.txt',header=FALSE)
+dim(xx)
+plot(xx[,3])
+plot(xx[,4])
+plot(xx[,3])
+plot(xx[1"970",3])
+plot(xx[1:970,3])
+write.table('foldbif.txt',xx[1:970,3])
+?write.table
+write.table(xx[1:970,3],'foldbif.txt')
+write.table(xx[1:970,3],'foldbif.txt',col.names=FALSE)
+write.table(xx[1:970,3],'foldbif.txt',col.names=FALSE,row.names=FALSE)
+xx=read.table('foldbif.txt',header=FALSE)
+plot(xx)
+xx
+dim(xx)
+plot(xx''
+plot(xx')
+cc=xx''
+cc=xx'
+plot(ts(xx))
+write.table
+?write.table
+xxx=data.frame(xx)
+dim(xxx)
+xxx
+head(xxx)
+xxx=data.frame(xx,col.names=FALSE)
+head(xxx)
+?data.frame
+xxx=data.frame(xx,row.names=NULL)
+head(xxx)
+write.table(xxx,'foldbif.txt',col.names=FALSE,row.names=FALSE)
+xx=read.table('foldbif.txt',header=FALSE)
+head(xx)
+xx=read.csv('foldbif.txt',header=FALSE)
+head(xx)
+?read.table
+xx=read.delim('foldbif.txt',header=FALSE)
+head(xx)
+plot(xx)
+xx=read.table('CSD_6Dec11.txt',header=FALSE)
+head(xx)
+plot(xx[,3])
+xx=read.delim('foldbif.txt',header=FALSE)
+head(xx)
+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")
+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
+?ch_ews
+?ch_ews
Added: pkg/R/bdstest_ews.R
===================================================================
--- pkg/R/bdstest_ews.R (rev 0)
+++ pkg/R/bdstest_ews.R 2012-01-17 16:12:59 UTC (rev 3)
@@ -0,0 +1,152 @@
+# BDS test Early Warning Signals
+# Author: Stephen R Carpenter, 22 Oct 2011
+# Modified by: Vasilis Dakos, January 1, 2012
+
+# BDS test bootstrapped on a time series
+BDSboot <- function(X,varname,nboot,epsvec,emb) { # begin function
+StdEpsAll <- X # name of variable for BDS
+neps <- length(epsvec)
+# Compute and print BDS test
+print('***********************************************',quote=FALSE)
+print(c('BDS test for ',varname),quote=FALSE)
+print(c('Embedding dimension = ',emb),quote=FALSE)
+BDS.data <- bds.test(StdEpsAll,m=emb,epsvec)
+print('BDS statistics for Nominal Data at each Epsilon',quote=FALSE)
+print(round(BDS.data$statistic,3))
+print('P value based on standard normal',quote=FALSE)
+print(round(BDS.data$p.value,3))
+# Bootstrap the BDS test
+nobs <- length(StdEpsAll)
+bootmat <- matrix(0,nrow=emb-1,ncol=neps) # matrix to count extreme BDS values
+for(i in 1:nboot) { # start bootstrap loop
+ epsboot <- sample(StdEpsAll,nobs,replace=TRUE)
+ BDS.boot <- bds.test(epsboot,m=emb,epsvec)
+ for(im in 1:(emb-1)) { # loop over embedding dimensions
+ bootvec <- BDS.boot$statistic[im,]
+ N.above <- ifelse(bootvec>BDS.data$statistic[im,],1,0)
+ bootmat[im,] <- bootmat[im,]+N.above
+ }
+ # Report progress: if hash is removed from the next two lines, the program
+ # will report each time an iteration is completed
+ # cat('iteration = ',i,' of ',nboot,'\n') #
+ # flush.console()
+ } # end bootstrap loop
+print(' ',quote=FALSE)
+print(c('Bootstrap P estimates for ',varname),quote=FALSE)
+print(c('Bootstrap iterations = ',nboot),quote=FALSE)
+Pboot <- bootmat/nboot
+for(im in 1:(emb-1)) {
+ print(c('For embedding dimension =',im+1),quote=FALSE)
+ print(c('For epsilon = ',round(epsvec,3),'bootstrap P = '),quote=FALSE)
+ print(Pboot[im,])
+ }
+print('**********************************************************',quote=FALSE)
+} # end function
+
+# MAIN FUNCTION
+bdstest_ews<-function(timeseries,ARMAoptim=TRUE,ARMAorder=c(1,0),embdim=3,epsilon=c(0.5,0.75,1),boots=1000,logtransform=FALSE,interpolate=FALSE){
+
+ # Load libraries
+ require('tseries')
+
+ 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)}
+
+ # Detrend the data
+ Eps1 <- diff(Y)
+
+ # Define BDS parameters
+ nboot <- boots
+ emb <- embdim # embedding dimension
+ eps.sd <- sd(Eps1)
+ epsvec <- round(eps.sd* epsilon,6)
+
+ # Run BDS with bootstrapping
+ BDSboot(Eps1,c('Detrended data'),nboot,epsvec,emb)
+
+ # Fit ARMA model based on AIC
+ if (ARMAoptim==TRUE){
+ arma=matrix(,4,5)
+ for (ij in 1:4){
+ for (jj in 0:4){
+ ARMA<-arima(Y, order = c(ij,0,jj),method="ML",include.mean = TRUE)
+ arma[ij,jj+1]=ARMA$aic
+ ind=which(arma==min(arma),arr.in=TRUE)
+ armafit<-arima(Y, order = c(ind[1],0,ind[2]-1),include.mean = TRUE)
+ print('ARMA model',quote=FALSE)
+ print(armafit,digits=4)
+ Eps2 <- armafit$residuals#[2:armafit$n.used]
+ }
+ }
+ }else{
+ armafit <- arima(Y, order = c(ARMAorder[1],0,ARMAorder[2]),include.mean = TRUE)
+ print('ARMA model',quote=FALSE)
+ print(armafit,digits=4)
+ Eps2 <- armafit$residuals#[2:armafit$n.used]
+ }
+
+ # Define BDS parameters
+ nboot <- boots
+ emb <- embdim # embedding dimension
+ eps.sd <- sd(Eps2)
+ epsvec <- round(eps.sd*epsilon,6)
+
+ # Run BDS with bootstrapping
+ BDSboot(Eps2,c('ARMA model residuals'),nboot,epsvec,emb)
+
+ # Fit GARCH(0,1) model to detrended data
+ Gfit <- garch(Y,order=c(0,1))
+ print('GARCH(0,1) model fit to detrended data',quote=FALSE)
+ print(Gfit,digits=4)
+
+ Eps3 <- Gfit$residuals[2:length(Y)]
+
+ # Define BDS parameters
+ nboot <- boots
+ emb <- embdim # embedding dimension
+ eps.sd <- sd(Eps3)
+ epsvec <- round(eps.sd*epsilon,6)
+
+ # Run BDS with bootstrapping
+ BDSboot(Eps3,c('GARCH(0,1) model residuals'),nboot,epsvec,emb)
+
+ # Plot the data
+ dev.new()
+ par(fig=c(0,0.5,0.5,1),mar=c(3, 4, 3, 2),cex.axis=0.8,cex.lab=1,mfrow=c(2,2),mgp=c(1.5,0.5,0),oma=c(1,1,2,1))
+ plot(timeindex,Y,type='l',col='black',lwd=1.7,xlab='time',ylab='data')
+ par(fig=c(0.5,1,0.8,0.95),mar=c(0, 4, 0, 2),new=TRUE)
+ plot(timeindex[1:(length(timeindex)-1)],Eps1,type="l",col="red",lwd=1,xlab="",ylab="",xaxt="n")
+ legend("topright","first-diff",bty="n",cex=0.8)
+ par(fig=c(0.5,1,0.65,0.8), new=TRUE)
+ plot(timeindex[1:(length(timeindex))],Eps2,type="l",xlab="",ylab="residuals",xaxt="n",col="green",lwd=1)
+ legend("topright","AR",bty="n",cex=0.8)
+ par(fig=c(0.5,1,0.5,0.65), new=TRUE)
+ plot(timeindex[1:(length(timeindex)-1)],Eps3,type="l",col="blue",xlab="time",ylab="",lwd=1)
+ legend("topright","GARCH",bty="n",cex=0.8)
+ mtext("time",side=1,outer=FALSE,line=1.4,cex=0.8)
+ # Time series diagnostics
+ par(fig=c(0,0.5,0,0.4),mar=c(3, 4, 0, 2), new=TRUE)
+ acf(Y,lag.max=25,main="")
+ par(fig=c(0.5,1,0,0.4), new=TRUE)
+ pacf(Y,lag.max=25,main="")
+ mtext("BDS_test Diagnostics",side=3,line=0.2, outer=TRUE)
+ }
\ No newline at end of file
Added: pkg/R/ch_ews.R
===================================================================
--- pkg/R/ch_ews.R (rev 0)
+++ pkg/R/ch_ews.R 2012-01-17 16:12:59 UTC (rev 3)
@@ -0,0 +1,88 @@
+# 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,logtransform=FALSE,interpolate=FALSE){
+ timeseries<-ts(timeseries) #strict data-types the input data as tseries object for use in later steps
+ if (dim(timeseries)[2]==1){
+ ts.in=timeseries
+ timeindex=1:dim(timeseries)[1]
+ }else if(dim(timeseries)[2]==2){
+ ts.in=timeseries[,2]
+ timeindex=timeseries[,1]
+ }else{
+ warning("not right format of timeseries input")
+ }
+
+ # Interpolation
+ if (interpolate){
+ YY<-approx(timeindex,ts.in,n=length(ts.in),method="linear")
+ ts.in<-YY$y
+ }
+
+ # Log-transformation
+ if (logtransform){
+ ts.in<-log(ts.in+1)}
+
+ 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)
+}
\ No newline at end of file
Added: pkg/R/ddjnonparam_ews.R
===================================================================
--- pkg/R/ddjnonparam_ews.R (rev 0)
+++ pkg/R/ddjnonparam_ews.R 2012-01-17 16:12:59 UTC (rev 3)
@@ -0,0 +1,153 @@
+# Drift Diffusion Jump Nonparametrics Early Warning Signals
+# Author: Stephen R Carpenter, 15 December 2011
+# Modified by: Vasilis Dakos, January 4, 2012
+
+# Function to compute Bandi, Johannes etc estimators for time series x
+# Inputs:
+# x0 is the regressor
+# dx is the first difference of x0
+# nx is number of first differences
+# DT is time step
+# bw is the bandwidth for the kernel
+# na is number of a values for computing the kernel
+# avec is the mesh for the kernel
+
+Bandi5 <- function(x0,dx,nx,DT,bw,na,avec) {
+# Set up constants and useful preliminaries
+SF <- 1/(bw*sqrt(2*pi)) # scale factor for kernel calculation
+x02 <- x0*x0 # second power of x
+dx2 <- dx*dx # second power of dx
+dx4 <- dx2*dx2 # fourth power of dx
+dx6 <- dx2*dx4 # sixth power of dx
+# Compute matrix of kernel values
+Kmat <- matrix(0,nrow=na,ncol=nx)
+for(i in 1:(nx)) { # loop over columns (x0 values)
+ Kmat[,i] <- SF*exp(-0.5*(x0[i]-avec)*(x0[i]-avec)/(bw*bw))
+ }
+# Compute M1, M2, M4, moment ratio and components of variance for each value of a
+M1.a <- rep(0,na)
+M2.a <- rep(0,na)
+M4.a <- rep(0,na)
+M6M4r <- rep(0,na) # vector to hold column kernel-weighted moment ratio
+mean.a <- rep(0,na) # centering of conditional variance
+SS.a <- rep(0,na) # sum of squares
+for(i in 1:na) { # loop over rows (a values)
+ Ksum <- sum(Kmat[i,]) # sum of weights
+ M1.a[i] <- (1/DT)*sum(Kmat[i,]*dx)/Ksum
+ M2.a[i] <- (1/DT)*sum(Kmat[i,]*dx2)/Ksum
+ M4.a[i] <- (1/DT)*sum(Kmat[i,]*dx4)/Ksum
+ M6.c <- (1/DT)*sum(Kmat[i,]*dx6)/Ksum
+ M6M4r[i] <- M6.c/M4.a[i]
+ mean.a[i] <- sum(Kmat[i,]*x0[2:(nx+1)])/Ksum
+ SS.a[i] <- sum(Kmat[i,]*x02[2:(nx+1)])/Ksum
+ }
+# Compute conditional variance
+S2.x <- SS.a - (mean.a*mean.a) # sum of squares minus squared mean
+# Compute jump frequency, diffusion and drift
+sigma2.Z <- mean(M6M4r)/(5) # average the column moment ratios
+lamda.Z <- M4.a/(3*sigma2.Z*sigma2.Z)
+sigma2.dx <- M2.a - (lamda.Z*sigma2.Z)
+# set negative diffusion estimates to zero
+diff.a <- ifelse(sigma2.dx>0,sigma2.dx,0)
+sigma2.dx <- M2.a # total variance of dx
+mu.a <- M1.a
+outlist <- list(mu.a,sigma2.dx,diff.a,sigma2.Z,lamda.Z,S2.x)
+# outputs of function:
+# mu.a is drift
+# sigma2.dx is total variance of dx
+# diff.a is diffusion
+# sigma2.Z is jump magnitude
+# lamda.Z is jump frequency
+# S2.x is conditional variance
+return(outlist)
+} # end Bandi function
+
+# MAIN FUNCTION
+ddjnonparam_ews<-function(timeseries,bandwidth=0.6,na=500,logtransform=TRUE,interpolate=FALSE){
+ require('KernSmooth')
+
+ 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)}
+
+ # Preliminaries
+ Xvec1<-Y
+ Tvec1<-timeindex
+ dXvec1 <- diff(Y)
+
+ DT <- Tvec1[2]-Tvec1[1]
+ bw <- bandwidth*sd(Xvec1) # bandwidth
+ alow <- min(Xvec1)
+ ahigh <- max(Xvec1)
+ na <- na
+ avec <- seq(alow,ahigh,length.out=na)
+ nx <- length(dXvec1)
+
+ # Bandi-type estimates
+ ParEst <- Bandi5(Xvec1,dXvec1,nx,DT,bw,na,avec)
+ Drift.vec <- ParEst[[1]]
+ TotVar.dx.vec <- ParEst[[2]]
+ Diff2.vec <- ParEst[[3]]
+ Sigma2Z <- ParEst[[4]]
+ LamdaZ.vec <- ParEst[[5]]
+ S2.vec <- ParEst[[6]]
+
+ # Interpolate time courses of indicators
+ TotVar.i <- approx(x=avec,y=TotVar.dx.vec,xout=Xvec1)
+ TotVar.t <- TotVar.i$y
+ Diff2.i <- approx(x=avec,y=Diff2.vec,xout=Xvec1)
+ Diff2.t <- Diff2.i$y
+ Lamda.i <- approx(x=avec,y=LamdaZ.vec,xout=Xvec1)
+ Lamda.t <- Lamda.i$y
+ S2.i <- approx(x=avec,y=S2.vec,xout=Xvec1)
+ S2.t <- S2.i$y
+
+ # Plot the data
+ dev.new()
+ par(mfrow=c(2,1),mar=c(3, 3, 2, 2),mgp=c(1.5,0.5,0),oma=c(1,1,1,1))
+ plot(Tvec1,Xvec1,type='l',col='black',lwd=2,xlab='',ylab='original data')
+ grid()
+ plot(Tvec1[1:length(Tvec1)-1],dXvec1,type='l',col='black',lwd=2,xlab='time',ylab='first-diff data')
+ grid()
+
+ # Plot indicators versus a
+ dev.new()
+ par(mfrow=c(2,2),mar=c(3, 3, 2, 2) ,cex.axis=1,cex.lab=1,mgp=c(2,1,0),oma=c(1,1,2,1))
+ plot(avec,S2.vec,type='l',lwd=1,col='black',xlab='a',ylab='conditional variance')
+ plot(avec,TotVar.dx.vec,type='l',lwd=1,col='blue',xlab='a',ylab='total variance of dx')
+ plot(avec,Diff2.vec,type='l',lwd=1,col='green',xlab='a',ylab='diffusion')
+ plot(avec,LamdaZ.vec,type='l',lwd=1,col='red',xlab='a',ylab='jump intensity')
+ mtext("DDJ nonparametrics versus a",side=3,line=0.1,outer=TRUE)
+
+ # Plot indicators versus time
+ dev.new()
+ par(mfrow=c(2,2),mar=c(3, 3, 2, 2),cex.axis=1,cex.lab=1,mgp=c(1.5,0.5,0),oma=c(1,1,2,1))
+ plot(Tvec1,S2.t,type='l',lwd=1,col='black',xlab='time',ylab='conditional variance')
+ plot(Tvec1,TotVar.t,type='l',lwd=1,col='blue',xlab='time',ylab='total variance of dx')
+ plot(Tvec1,Diff2.t,type='l',lwd=1,col='green',xlab='time',ylab='diffusion')
+ plot(Tvec1,Lamda.t,type='l',lwd=1,col='red',xlab='time',ylab='jump intensity')
+ mtext("DDJ nonparametrics versus time",side=3,line=0.1,outer=TRUE)
+
+ # Output
+ nonpar_x<-data.frame(avec,S2.vec,TotVar.dx.vec,Diff2.vec,LamdaZ.vec)
+ nonpar_t<-data.frame(Tvec1,S2.t,TotVar.t,Diff2.t,Lamda.t)
+ return(c(nonpar_x,nonpar_t))
+ }
\ No newline at end of file
Added: pkg/R/earlywarnings-internal.R
===================================================================
--- pkg/R/earlywarnings-internal.R (rev 0)
+++ pkg/R/earlywarnings-internal.R 2012-01-17 16:12:59 UTC (rev 3)
@@ -0,0 +1,126 @@
+.Random.seed <-
+c(403L, 10L, 1897936630L, -1719655904L, 881993699L, -1426011007L,
+-545453328L, 1718534542L, -512610567L, 806804619L, 1551647930L,
+1053810620L, 1911463167L, -1908805547L, 1750049644L, 2123865026L,
+-142685155L, -1759859561L, -1014494722L, -1217082952L, 1255720347L,
+2006477817L, -886901944L, -309207434L, -460957775L, 745187811L,
+331939506L, 1853436676L, -2002112441L, 533583645L, 889597332L,
+642175290L, 1125657989L, 604317279L, 400494822L, 178410352L,
+-82086765L, 1122496817L, -2036508704L, -1899635586L, -1732489463L,
+1004474683L, 331183370L, 1501129964L, 135584111L, 553475749L,
+1369354332L, 4532274L, -876093651L, 1237888487L, -689442546L,
+-1154345240L, -522045205L, 315701641L, 1326041368L, -903521210L,
+679208097L, 478710579L, -1902939486L, 913820052L, -885991401L,
+766689293L, 1430324516L, 1340555530L, -2125858699L, -1786883121L,
+-253631274L, -1929153408L, -1662547261L, -166800287L, 1767019728L,
+339002926L, 1357395929L, 242964459L, 739249754L, 2142111772L,
+873600287L, -1285188235L, 1939907596L, 1875764706L, -1895803203L,
+1101754679L, 2056639966L, -491255016L, -762183045L, 317141017L,
+-1126897496L, 588696150L, 769111761L, -1558038909L, -1884820206L,
+-996247388L, 318617703L, -500631107L, -1152445644L, 1834481818L,
+744742629L, 1909700735L, 1331229190L, -842179184L, 1398014003L,
+-11243567L, -1610422912L, -1080301474L, 1667155881L, -20775333L,
+1891000938L, -162574580L, -621156017L, -480783867L, 2044100796L,
+-1994020462L, -360136051L, -1822199609L, -1733360210L, 741478408L,
+-1374419957L, 469090537L, 1340784568L, -1498539930L, -145405695L,
+-667746413L, 298291266L, -160737292L, 1008992759L, 646776301L,
+1482966660L, -2073347926L, -754481643L, -1261741393L, -265678666L,
+-722563360L, -2146725341L, -1214038719L, -1085439952L, -66938290L,
+-580264263L, -408507701L, -2041156742L, -396322180L, -1975112001L,
+-1247111275L, 649809836L, 716475778L, 413700957L, 1484275287L,
+63397182L, 1189424504L, 2144091739L, 1594251705L, -123671544L,
+-705005898L, 13796849L, -810116061L, 47248498L, 1563797188L,
+1056604295L, 2101675101L, -2003195308L, -2060451718L, 624264133L,
+878130463L, 1143837862L, 1336550960L, 101940435L, -1417497871L,
+189612576L, -2088655682L, 142430409L, -461516165L, -1500874678L,
+-860059476L, -107614289L, 1220885349L, 1778667804L, 246208882L,
+1748781549L, 2065243687L, -1784066482L, 1987801768L, 1644318891L,
+-1741915703L, 1497126360L, 359048326L, 896649569L, -1930998541L,
+-1583516830L, -1234928172L, -847497001L, 995543757L, -1258546588L,
+-1743007670L, -103279179L, -2136114929L, 1909679638L, -2092556096L,
+-645581181L, -1269108063L, -1472605040L, 993591150L, -1727801063L,
+-2075838037L, 1775802266L, 112658268L, 1742345823L, -681433547L,
+-1125923380L, -580319966L, -1233211523L, 333218935L, -1102366306L,
+1673211992L, -1803375173L, -633584807L, 334377192L, -2146692842L,
+835498641L, 1192542019L, -31352238L, 576433380L, -1146865881L,
+-605998211L, -2042247564L, 1718532186L, -1210171227L, -736632129L,
+1835678022L, 1496747984L, -1934177805L, -1679943776L, -354409284L,
+-414342704L, -133320350L, 1767960552L, 961975244L, 385066748L,
+1152828402L, -532697984L, 2003863796L, -1274437496L, -1399242438L,
+-1537692848L, 20545644L, 536116372L, -1469086382L, 2000079584L,
+-439759668L, -204983136L, 1194440066L, 1692673256L, -1416100628L,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/earlywarnings -r 3
More information about the Earlywarnings-commits
mailing list