rollWarn: Roll Warn Calculate early warnings statistics over a rolling...

In rBatt/fishyWarnings: How should management respond to warning signals of haddock collapse?

Description

Roll Warn Calculate early warnings statistics over a rolling window

Usage

rollWarn(x, stat = eval(formals(ews)$stat), win = 15)

Arguments

x	numeric vector over which to calculate statistics; should be a regularly spaced time series
stat	name of statistic to use; see ews
win	window length

Value

numeric vector of a rolling window statistic

See Also

ews

Examples

# ==========================================
# = Set up simulation and analysis options =
# ==========================================
nYear <- 2E2
Year <- 1:nYear
tripVal <- c(0.36, 0.21)
nBurn <- 50
sdU <- 0.05

# set up qE to increase harvest until collapse
qE_range <- c(tripVal[1]-0.2, tripVal[1]+0.05)

# goal here is to have harvest rate precisely hit certain
# values of qE during an integer year,
# while retaining linear change throughout
qEPoints <- c(qE_range[1], tripVal[2], tripVal[1], qE_range[2])
qEPoints <- qEPoints[qEPoints>=min(qE_range) & qEPoints<=max(qE_range)]
seq_arg0 <- c(as.list(range(qEPoints)), list(length.out=nYear))
yrProbs <- (qEPoints-min(qEPoints))/diff(range(qEPoints))
yrPoints <- quantile(Year, yrProbs)
qE <- approx(x=yrPoints, y=qEPoints, xout=Year)$y

# ==============================
# = Claculate time series of B =
# ==============================
B0 <- 400
B0 <- Burn(B0, n=nBurn, qE=qE[1], sdU=0)
Bvec <- c(B0, rep(NA, nYear-1))
for(i in 2:nYear){
	Bvec[i] <- Bstep(B=Bvec[i-1], qE=qE[i], sdU=sdU)
}

# ==================================
# = Calculate rolling window stats =
# ==================================
sdVec <- rollWarn(Bvec, stat='sd', win=min(50, nYear/5))
ac1Vec <- rollWarn(Bvec, stat='ac1', win=min(50, nYear/5))
redVec0 <- rollWarn(Bvec, stat="redShift", win=min(50, nYear/5)) # holds the full spectrum
redList <- attributes(redVec0)$rsl # this gets the full spectrum formatted for plotting
redVec <- sapply(redVec0, function(x)x[[1]]) # this gets just the spectral slope

# =======================
# = Plot B and Warnings =
# =======================
par(mfrow=c(2,2), mar=c(1.75,1.75,0.5,1.75), ps=8, cex=1, mgp=c(0.75,0.15,0), tcl=-0.15)
plot(Year, Bvec, type='l')
abline(v=yrPoints, lty='dashed')
text_y <- min(Bvec)+yrProbs*diff(range(Bvec))
text(yrPoints, y=text_y, label=paste0("qE=", qEPoints), pos=c(4,4,2,2))

plot(Year, sdVec, col="forestgreen", type="l", ylab="Standard Deviation")
abline(v=yrPoints, lty='dashed')
plot(Year,ac1Vec, col="blue", type='l', ylab="AR(1)")
abline(v=yrPoints, lty='dashed')
plot(redList, xaxs='r')
par(new=TRUE)
plot(Year, redVec, col=adjustcolor("white",0.5), lwd=3, type='l', ylab="", xaxt='n', yaxt='n')
lines(Year, redVec, col=adjustcolor("black",0.5), lwd=0.5)
axis(side=4)
mtext("redShift", side=4, line=0.75)
abline(v=yrPoints, lty='dashed')

rBatt/fishyWarnings documentation built on May 26, 2019, 7:45 p.m.