Skip to content
GitLab
Commit d77076fc authored by Iago Mosqueira's avatar Iago Mosqueira
Browse files

Cleaning

parent 1d12514c
Expand all Hide whitespace changes
Inline Side-by-side
DATA.bib 0 → 100644
View file @ d77076fc
@Misc{dinum.txt,
originator = {WGNSSK},
year = {2020},
title = {Discards by age},
period = {1957-2018},
access = {Public},
source = {file},
}
@Misc{fleet.txt,
originator = {WGNSSK},
year = {2020},
title = {BTS and SNS survey indices by age},
period = {},
access = {Public},
source = {file},
}
@Misc{fprop.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{index.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{lanum.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{laton.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{matprop.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{mprop.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{natmor.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{new,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{wedi.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{wela.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
@Misc{west.txt,
originator = {WGNSSK},
year = {2020},
title = {},
period = {},
access = {Public},
source = {file},
}
SOFTWARE.bib 0 → 100644
View file @ d77076fc
@Manual{AAP,
author = {Iago Mosqueira},
year = {2020},
title = {AAP: Aarts and Poos Stock Assessment Model that Estimates Bycatch},
version = {0.1.1, branch feature/knot_branch, last commit 5 Feb 2020},
source = {iagomosqueira/AAP@e79077e},
}
TMP_SAM.R 0 → 100644
View file @ d77076fc
# modelSAM.R - DESC
# sol.27.4-sa/modelSAM.R
# Copyright Niels HINTZEN (WMR), 2020
# Author: Niels HINTZEN (WMR) <niels.hintzen@wur.nl>
# Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
#
# Distributed under the terms of the EUPL-1.2
# devtools::install_github('fishfollower/SAM/stockassessment', ref='components')
library(FLSAM)
load('data/data.RData')
#- Correct some missing data in SOL data
SOL <- stock
# SOL@catch.wt[ac(10),ac(2019)] <- SOL@landings.wt[ac(10),ac(2018)]
# SOL@discards.n[ac(10),ac(2019)] <- 0
# SOL@discards.wt[ac(10),ac(2019)] <- 0
# SOL@stock.wt[,ac(2019)] <- rowMeans(SOL@stock.wt[,ac(2016:2018)])
# SOL@m[,ac(2019)] <- SOL@m[,ac(2018)]
# SOL@mat[,ac(2019)] <- SOL@mat[,ac(2018)]
# SOL@m.spwn[,ac(2019)] <- SOL@m.spwn[,ac(2018)]
# SOL@harvest.spwn[,ac(2019)] <- SOL@harvest.spwn[,ac(2018)]
#
#- Correct some settings in tuning data
SOL.tun <- indices[c("GAM", "SNS")]
type(SOL.tun[[1]]) <- "number"
type(SOL.tun[[2]]) <- "number"
#- Setup ctrl file
SOL.ctrl <- FLSAM.control(SOL,SOL.tun)
SOL.ctrl@residuals <- FALSE
#- Run default assessment
SOL.sam <- FLSAM(SOL,SOL.tun,SOL.ctrl)
defAIC <- AIC(SOL.sam); defAIC
#---------------------
#- Add correlation structure for surveys
#---------------------
SOL.ctrl@cor.obs[2,1:8] <- c(101:107,107)
SOL.ctrl@cor.obs[3,1:5] <- c(201:204,204)
SOL.ctrl@cor.obs.Flag[2:3] <- as.factor("AR")
SOL.ctrl <- update(SOL.ctrl)
SOL.samcorObs <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.sam)
corObsAIC <- AIC(SOL.samcorObs); corObsAIC
xyplot(value + lbnd + ubnd ~ age | fleet,data=cor.obs(SOL.samcorObs),type="l",col=c(1,"grey","grey"),lty=c(2,1,1),ylim=c(0,5))
#- Decision
SOL.ctrl@cor.obs[2,1:8] <- c(rep(101,2),rep(102,6))
SOL.ctrl@cor.obs[3,1:5] <- c(rep(201,2),rep(202,3))
SOL.samcorObs <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.sam)
corObsAIC <- AIC(SOL.samcorObs); corObsAIC
#---------------------
#- Configure F random walks
#---------------------
SOL.ctrl@f.vars["catch unique",] <- c(1:9,9)
SOL.ctrl <- update(SOL.ctrl)
SOL.samfvar <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samcorObs)
# DEBUG Error in solve.default(h, g) :
# Lapack routine dgesv: system is exactly singular: U[20,20] = 0
fvarAIC <- AIC(SOL.samfvar); fvarAIC
xyplot(value + lbnd + ubnd ~ age | fleet,data=f.var(SOL.samfvar),type="l",col=c(1,"grey","grey"),lty=c(2,1,1))
#- Decision
SOL.ctrl@f.vars["catch unique",] <- c(1,2,rep(3,6),rep(4,2))
SOL.ctrl <- update(SOL.ctrl)
SOL.samfvar <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samcorObs)
fvarAIC <- AIC(SOL.samfvar); fvarAIC
#---------------------
#- Configure F correlation structure (usually sensible to also run retrospectives)
#---------------------
SOL.ctrl@cor.F <- 2
SOL.samcorF2 <- SOL.samfvar
SOL.retrocorF2 <- retro(SOL,SOL.tun,SOL.ctrl,retro=5)
SOL.ctrl@cor.F <- 1
SOL.samcorF1 <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samfvar)
SOL.retrocorF1 <- retro(SOL,SOL.tun,SOL.ctrl,retro=5)
SOL.ctrl@cor.F <- 0
SOL.samcorF0 <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samfvar)
SOL.retrocorF0 <- retro(SOL,SOL.tun,SOL.ctrl,retro=5)
corFAICs <- AIC(FLSAMs(corF0=SOL.samcorF0,corF1=SOL.samcorF1,corF2=SOL.samcorF2)); corFAICs
storeMohnsRho <- matrix(NA,nrow=3,ncol=3,dimnames=list(type=c("ssb","fbar","rec"),model=c("corF2","corF1","corF0")))
storeMohnsRho["ssb","corF2"] <- mean(mohns.rho(SOL.retrocorF2,type="ssb",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["fbar","corF2"] <- mean(mohns.rho(SOL.retrocorF2,type="fbar",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["rec","corF2"] <- mean(mohns.rho(SOL.retrocorF2,type="rec",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["ssb","corF1"] <- mean(mohns.rho(SOL.retrocorF1,type="ssb",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["fbar","corF1"] <- mean(mohns.rho(SOL.retrocorF1,type="fbar",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["rec","corF1"] <- mean(mohns.rho(SOL.retrocorF1,type="rec",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["ssb","corF0"] <- mean(mohns.rho(SOL.retrocorF0,type="ssb",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["fbar","corF0"] <- mean(mohns.rho(SOL.retrocorF0,type="fbar",ref.year=2019,span=5)[1:5,1])
storeMohnsRho["rec","corF0"] <- mean(mohns.rho(SOL.retrocorF0,type="rec",ref.year=2019,span=5)[1:5,1])
storeMohnsRho
#- Decision
SOL.ctrl@cor.F <- 2
#---------------------
#- Configure Observation variances
#---------------------
SOL.ctrl@obs.vars["catch unique",] <- c(1:9,9)
SOL.ctrl@obs.vars["BTS-ISIS",ac(1:9)] <- c(101:108,108)
SOL.ctrl@obs.vars["SNS",ac(1:6)] <- c(201:205,205)
SOL.ctrl <- update(SOL.ctrl)
SOL.samobsVar <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samcorF2)
obsVarAIC <- AIC(SOL.samobsVar); obsVarAIC
xyplot(value + lbnd + ubnd ~ age | fleet,data=obs.var(SOL.samobsVar),type="l",col=c(1,"grey","grey"),lty=c(2,1,1))
#- Decision
SOL.ctrl@obs.vars["catch unique",] <- c(1,2,3,3,rep(4,6))
SOL.ctrl@obs.vars["BTS-ISIS",ac(1:9)] <- c(rep(101,5),rep(102,2),rep(103,2))
SOL.ctrl@obs.vars["SNS",ac(1:6)] <- c(rep(201,3),rep(202,3))
SOL.ctrl <- update(SOL.ctrl)
SOL.samobsVar <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samcorF2)
obsVarAIC <- AIC(SOL.samobsVar); obsVarAIC
xyplot(value + lbnd + ubnd ~ age | fleet,data=obs.var(SOL.samobsVar),type="l",col=c(1,"grey","grey"),lty=c(2,1,1))
#---------------------
#- Configure catchabilities
#---------------------
xyplot(value + lbnd + ubnd ~ age | fleet,data=catchabilities(SOL.samobsVar),type="l",col=c(1,"grey","grey"),lty=c(2,1,1),scales=list(y="free"))
#- Decision
SOL.ctrl@catchabilities["BTS-ISIS",ac(1:9)] <- c(1,1,2,3,4,5,5,6,6)
SOL.ctrl@catchabilities["SNS",ac(1:6)] <- c(1,2,3,4,4,4) + 101
SOL.ctrl <- update(SOL.ctrl)
SOL.samcatch <- FLSAM(SOL,SOL.tun,SOL.ctrl,starting.values=SOL.samobsVar)
catchAIC <- AIC(SOL.samcatch); catchAIC
#-------------------------------------------------------------------------------
#- Run final model
#-------------------------------------------------------------------------------
SOL.ctrl@residuals <- TRUE
SOL.sam <- FLSAM(SOL,SOL.tun,SOL.ctrl)
SOL.ctrl@residuals <- FALSE
SOL.retro <- retro(SOL,SOL.tun,SOL.ctrl,retro=5)
SOL.looi <- looi(SOL,SOL.tun,SOL.ctrl,type="full")
save(SOL,SOL.tun,SOL.ctrl,SOL.sam,SOL.retro,SOL.looi,
file=file.path(path, "sam.RData"))
#- Alternative to fit directly into stockassessment
data <- FLSAM2SAM(FLStocks(residual=stock),tune)
conf <- ctrl2conf(control, data)
par <- stockassessment::defpar(data,conf)
fit <- stockassessment::sam.fit(data,conf,par)
rsd <- residuals(fit)
#-------------------------------------------------------------------------------
# Do the plotting
#-------------------------------------------------------------------------------
pdf(file.path("model", "sam", "plots_diagnostics.pdf"))
print(plot(SOL.sam,futureYrs=F))
residual.diagnostics(SOL.sam)
resids <- residuals(SOL.sam)
resids$std.res[which(is.na(resids$std.res))] <- 0
print(xyplot(age ~ year | fleet,data=resids,main="Residuals by fleet",group=resids$fleet,cex=resids$std.res,
panel=function(...){
lst <- list(...)
panel.xyplot(lst$x,lst$y,pch=ifelse(lst$cex[lst$subscript]>0,1,19),col="black",cex=1*abs(lst$cex[lst$subscript]))
}))
print(xyplot(age ~ fleet | as.factor(year),data=resids,main="Residuals by year",group=resids$fleet,cex=resids$std.res,scales=list(x=list(rot=90)),
panel=function(...){
lst <- list(...)
panel.xyplot(lst$x,lst$y,pch=ifelse(lst$cex[lst$subscript]>0,1,19),col="black",cex=1*abs(lst$cex[lst$subscript]))
}))
# figure - catchabilities at age from HERAS
catch <- catchabilities(SOL.sam)
print(xyplot(value+ubnd+lbnd ~ age | fleet,catch,
scale=list(alternating=FALSE,y=list(relation="free")),as.table=TRUE,
type="l",lwd=c(2,1,1),col=c("black","grey","grey"),
subset=fleet %in% names(SOL.tun),
main="Survey catchability parameters",ylab="Catchability",xlab="Age"))
# figure - variance by data source
obv <- obs.var(SOL.sam)
obv$str <- paste(obv$fleet,ifelse(is.na(obv$age),"",obv$age))
obv <- obv[order(obv$value),]
bp <- barplot(obv$value,ylab="Observation Variance",
main="Observation variances by data source",col=factor(obv$fleet))
axis(1,at=bp,labels=obv$str,las=3,lty=0,mgp=c(0,0,0))
legend("topleft",levels(obv$fleet),pch=15,col=1:nlevels(obv$fleet),pt.cex=1.5)
# figure - variance vs uncertainty for each data source
plot(obv$value,obv$CV,xlab="Observation variance",ylab="CV of estimate",log="x",
pch=16,col=obv$fleet,main="Observation variance vs uncertainty")
text(obv$value,obv$CV,obv$str,pos=4,cex=0.75,xpd=NA)
# figure - fishing age selectivity per year
sel.pat <- merge(f(SOL.sam),fbar(SOL.sam),
by="year",suffixes=c(".f",".fbar"))
sel.pat$sel <- sel.pat$value.f/sel.pat$value.fbar
sel.pat$age <- as.numeric(as.character(sel.pat$age))
print(xyplot(sel ~ age|sprintf("%i's",floor((year)/5)*5),sel.pat,
groups=year,type="l",as.table=TRUE,
scale=list(alternating=FALSE),
main="Selectivity of the Fishery by Pentad",xlab="Age",ylab="F/Fbar"))
# figure - correlation matrix of model parameters
print(cor.plot(SOL.sam))
#Plot uncertainties as a function of time
CV.yrs <- ssb(SOL.sam)$year
CV.dat <- cbind(SSB=ssb(SOL.sam)$CV,
Fbar=fbar(SOL.sam)$CV,Rec=rec(SOL.sam)$CV)
print(matplot(CV.yrs,CV.dat,type="l",ylim=range(pretty(c(0,CV.dat))),yaxs="i",
xlab="Year",ylab="CV of estimate",main="Uncertainties of key parameters"))
legend("topleft",legend=colnames(CV.dat),lty=1:5,col=1:6,bty="n")
print(plot(SOL.retro))
retroParams(SOL.retro)
yrs <- 2010:2019
res <- lapply(SOL.retro, f)
res <- lapply(res, function(y) {
y[which(y$year %in% (max(y$year) - 20):(max(y$year))), ]
})
res <- lapply(res, function(y) {
cbind(y, retro = max(y$year))
})
res <- do.call(rbind, res)
res <- subset(res, year %in% yrs)
print(xyplot(value ~ an(age) | as.factor(year), data = res,
type = "l", groups = retro, auto.key = list(space = "right",
points = FALSE, lines = TRUE, type = "l"), main = paste("Retrospective pattern in F at age"),
ylab = "F", xlab = "Ages", panel = panel.superpose, panel.groups = function(...) {
panel.grid(v = -1, h = -1, lty = 3)
panel.xyplot(...)
}, scales = list(alternating = 1, y = list(relation = "free",
rot = 0))))
print(plot(SOL.looi,main="Leave one in"))
dev.off()
bootstrap/initial/sa/aap2019.R 0 → 100644
View file @ d77076fc
# sa.R - DESC
# /sa.R
# Copyright Iago MOSQUEIRA (WMR), 2020
# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
#
# Distributed under the terms of the EUPL-1.2
library(AAP)
# --- WGNSSK 2019
wg <- mget(load("bootstrap/initial/sa/aap_wgnssk2019.RData"))
aap2019 <- stock
harvest(aap2019) <- wg$res$harvest
units(harvest(aap2019)) <- "f"
stock.n(aap2019) <- wg$res$stock.n
save(aap2019, file="aap2019.RData")
bootstrap/survey/BTS-GAM-age10plus.csv
View file @ d77076fc
This diff is collapsed.
data.R
View file @ d77076fc
......@@ -4,7 +4,7 @@
# Copyright Iago MOSQUEIRA (WMR), 2019
# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
#
# Distributed under the terms of the GPL 3.0
# Distributed under the terms of the EUPL-1.2
# INPUT: bootstrap/initial/data (VPA format)
# OUPUT: data/data.RData (FLR), data/*.csv (TAF tables)
......@@ -126,3 +126,4 @@ indextables <- lapply(FLQuants(
survey.isis=index(indices[["GAM-ISIS"]])), function(x) plus(flr2taf(x)))
write.taf(c(stocktables, stockfulltables, indextables), dir="data")
explore.R
View file @ d77076fc
......@@ -6,6 +6,33 @@
#
# Distributed under the terms of the EUPL-1.2
# --- a4a submodels
# fmodel
fmod <- ~te(age, year, k = c(7, 22), bs = "tp")
# cohort
fmod <- ~s(replace(age, age > 8, 8), k = 4) +
s(pmax(year - age, 1957), k = 8) + s(year, k = 22)
# AAP-like
fmod <- ~te(replace(age, age > 8, 8), year, k = c(7, 22), bs = "tp")
fmod <- ~s(replace(age, age > 8, 8), k=6) +
s(year, k=22, by=breakpts(age, c(2, 6, 8)))
# srmodel
srmod <- ~factor(year)
# vmodel (catch, GAM, SNS)
# SAM-like
vmod <- list(
~s(age, k=4),
~breakpts(age, c(2, 5, 8)),
~breakpts(age, 3)
)
# --- PG runs
......@@ -100,76 +127,22 @@ dev.off()
run <- runs[["BTS-GAM"]]
stock <- stocks[["BTS-GAM"]]
ci <- function(name, std) {
x <- std[std$name==name,]
data.frame(lo=x$mean - 2 * x$stddev, hi=x$mean + 2 * x$stddev)
}
err <- data.table(
Run="gam",
year=1957:2018,
Rec=c(rec(stock)),
Rec_lo=exp(ci("log_initpop", run@stdfile)[-(1:9),])$lo,
Rec_hi=exp(ci("log_initpop", run@stdfile)[-(1:9),])$hi,
SSB=c(ssb(stock)),
SSB_lo=ci("SSBo", run@stdfile)$lo,
SSB_hi=ci("SSBo", run@stdfile)$hi,
Fbar=c(fbar(stock)),
Fbar_lo=ci("Fbar", run@stdfile)$lo,
Fbar_hi=ci("Fbar", run@stdfile)$hi)
library(patchwork)
p1 <- ggplot(err, aes(x=year)) +
geom_line(aes(y=SSB)) + ylim(c(0,NA)) +
geom_ribbon(aes(ymin=SSB_lo, ymax=SSB_hi), alpha=0.3, colour="grey")
mets <- metricsAAP(run)
p2 <- ggplot(err, aes(x=year)) +
geom_line(aes(y=Fbar)) + ylim(c(0,NA)) +
geom_ribbon(aes(ymin=Fbar_lo, ymax=Fbar_hi), alpha=0.3, colour="grey")
p3 <- ggplot(err, aes(x=year)) +
geom_line(aes(y=Rec)) + ylim(c(0,NA)) +
geom_ribbon(aes(ymin=Rec_lo, ymax=Rec_hi), alpha=0.3, colour="grey")
retr <- cbind(rbindlist(lapply(retro, function(x) data.table(model.frame(metrics(x,
list(SSB=ssb, Fbar=fbar, Rec=rec)), drop=TRUE))), idcol="Run"),
Rec_lo=NA, Rec_hi=NA, SSB_lo=NA, SSB_hi=NA, Fbar_lo=NA, Fbar_hi=NA)
retrerr <- rbind(err, retr)
rp1 <- ggplot(retrerr, aes(x=year)) +
geom_line(aes(y=SSB, group=Run, colour=Run)) +
geom_ribbon(data=retrerr[Run=="gam"], aes(ymin=SSB_lo, ymax=SSB_hi),
alpha=0.3, colour="grey") +
theme(legend.position="none") + ylim(c(0,NA))
rp2 <- ggplot(retrerr, aes(x=year)) +
geom_line(aes(y=Fbar, group=Run, colour=Run)) +
geom_ribbon(data=retrerr[Run=="gam"], aes(ymin=Fbar_lo, ymax=Fbar_hi),
alpha=0.3, colour="grey") +
theme(legend.position="none") + ylim(c(0,NA))
rp3 <- ggplot(retrerr, aes(x=year)) +
geom_line(aes(y=Rec, group=Run, colour=Run)) +
geom_ribbon(data=retrerr[Run=="gam"], aes(ymin=Rec_lo, ymax=Rec_hi),
alpha=0.3, colour="grey") +
theme(legend.position="none") + ylim(c(0,NA))
pdf(file="model/aap/model_aap.pdf")
p1 / p2 / p3
plot(runs$McMC + stock, metrics=list(SSB=ssb, Fbar=fbar, Rec=rec))
rp1 / rp2 / rp3
dev.off()
colnames(mets)[3] <- "data"
rmets <- rbindlist(lapply(retro, function(x) data.table(as.data.frame(metrics(x,
list(SSB=ssb, F=fbar, Rec=rec)), drop=TRUE))), idcol="run")
dat <- rbind(cbind(mets, run="base"), cbind(rmets, lowq=NA, uppq=NA))
ggplot(dat, aes(x=year)) +
geom_line(aes(y=data, colour=run, group=run)) + ylim(c(0,NA)) +
geom_ribbon(data=dat[run == "base", ],
aes(ymin=lowq, ymax=uppq), alpha=0.3, colour="grey") +
facet_grid(qname~., scales="free_y")
# LANDINGS.N
# --- LANDINGS.N
x <- catch.n(stock)[, ac(2011:2018)]
......
forecast.R 0 → 100644
View file @ d77076fc
# forecast.R - DESC
# /forecast.R
# Copyright Iago MOSQUEIRA (WMR), 2020
# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
#
# Distributed under the terms of the EUPL-1.2
# --- RUN short term forecast
# REC = geomean()
# REC = rct3()
model.R
View file @ d77076fc
# model.R - DESC
# /model.R
# Copyright Iago MOSQUEIRA (WMR), 2019
# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
# Distributed under the terms of the GPL 3.0
# Distributed under the terms of the EUPL-1.2
library(AAP)
library(icesAdvice)
library(icesTAF)
mkdir("model/aap")
load('data/data.RData')
# RUN models
set.seed(1844)
# --- RUN model
control <- AAP.control(pGrp=TRUE, qplat.surveys=8, qplat.Fmatrix=9,
Sage.knots=6, Fage.knots=8, Ftime.knots=28, Wtime.knots=5, mcmc=FALSE)
base <- aap(stock, indices[c("BTS-ISIS", "SNS")], control=control)
# GAM10p + SNS
fit <- aap(stock, indices[c("GAM10p", "SNS")], control=control,
stdfile=base@stdfile)
# --- RUN McMC
mcmc <- aap(stock, indices[c("GAM10p", "SNS")],