# Load data 
rm(list = ls())
Catch <- read.csv('Catch.csv')
cpue <- read.csv('cpue.csv')

myK <- 355.4
myr <- 0.104

# Setup data 
datstan <- list(N = nrow(Catch),
                LnCPUE = log(cpue$CPUEi),
                Catch = Catch$Catch,
                LnK_mu = log(myK),
                LnK_sd = 0.18,
                phi = 0.68,
                Lnr_mu = log(myr),
                Lnr_sd = 0.22
                )
inits <- list(
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.99, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.95, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.90, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.85, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.80, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.75, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.70, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*1.1), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*1.2), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*1.3), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*1.4), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*0.9), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*0.8), Lnq1 = log(1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*0.8), Lnq1 = log(0.9), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912),
  list(Lnr = log(myr), LnK = log(myK*0.8), Lnq1 = log(1.1), sigma = 0.2, tau = 0.1, phi = 0.65, m = 1.912)
)
# Run
library(rstan)
options(mc.cores = 16)
fitm1 <- stan(file = 'BSPM.stan', data = datstan,
              iter=40000, chains=16, thin = 40,
              control = list(max_treedepth = 24, stepsize = 0.5, 
                             stepsize_jitter = 0.1), 
              init = inits)