#########################################################################
# Parameter uncertainty attempts #
#########################################################################
multid_SDP_matrix <- function(f, p, x_grid, h_grid, sigma_g){
gridsizes <- c(a=5, s=6, x=7 ,h=4)
sigma_g <- 0.2
a_grid <- seq(1, 2, length=gridsizes["a"])
s_grid <- seq(.0, .1, length=gridsizes["s"])
x_grid <- seq(0, 1, length=gridsizes["x"])
h_grid <- seq(0, .5, length=gridsizes["h"])
f <- function(x,a,h,p){
max(0, (x - h) * a / (1 + (x - h) * p[1]))
}
p = c(1)
## For each value of x_i, find the probability of going to x_j (given a_i and s_i)
# P( x_i -> x_j; | a_i, s_i)
px <- function(a_i, s_i, x_i)
sapply(x_grid, function(x_j){
sum(
sapply(a_grid, function(a){
dlnorm(x_j, log(f(x_i, a, h, p)) - sigma_g ^ 2 / 2, sigma_g) *
dnorm(a, a_i, s_i)
})
)
})
x_
px(a_i, s_i, x_t_minus_1)
# a_i -> a_j | x_i, s_i
P_to_aj_given <-
sapply(a_grid, function(a_j){
dlnorm( log(f(x_tminus1_i, a_j, h,p) - sigma_g ^ 2 / 2, sigma_g)) * dnorm(a_j, a_i, s_i)
})
P_to_sj_given <-
dat <-
sapply(h_grid, function(h){
sapply(a_grid, function(a_hat){
sapply(s_grid, function(sigma_a){
sapply(x_grid, function(x){
sapply(a_grid, function(a){
sapply(x_grid, function(y){
dlnorm(y, log(f(x, a, h, p))-sigma_g^2/2, sigma_g) *
dnorm(a, a_hat, sigma_a)
})
})
})
})
})
})
get_column <- function(a,s,x,h, dat, gridsizes){
stride_a <- gridsize["x"] * gridsizes["a"] * gridsizes["x"] * gridsizes["s"]
stride_s <- gridsize["x"] * gridsizes["a"] * gridsizes["x"]
stride_x <- gridsize["x"] * gridsizes["a"]
startpt <-
1 + (a-1) * stride_a +
1 + (s - 1) * stride_s +
1 + (x - 1) * stride_x
endpt <- startpt + gridsizes["a"]
dat[startpt:endpt, h]
}
get_SDP_matrix <- function(){
n <- dim(dat)[1]
SDP_Mat <- matrix(NA, ncol=n, nrow=n)
a1s1x1
a1s1x1 get_column(1,1,1, 1, dat, gridsizes)
}
# Format of h_a_s_x_A:
# h1, h2 h3
# a1,x1,s1,a1
# a2,x1,s1
# a3,x1,s1
# a1,x2,s1
# a2,x2,s1
# a3 x2,s1
# a1 x3,s1
# ...
# a1,x1,s2
find_dp_optim <- function(SDP_Mat, x_grid, h_grid, OptTime, xT, profit,
delta, reward=0){
## Initialize space for the matrices
V <- rep(0,gridsize) # initialize BC,
# give a fixed reward for having value larger than xT at the end.
V[x_grid >= xT] <- reward # a "scrap value" for x(T) >= xT
# loop through time
for(time in 1:OptTime){
# try all potential havest rates
V1 <- sapply(1:HL, function(i){
# Transition matrix times V gives dist in next time
SDP_Mat[[i]] %*% V +
# then (add) harvested amount times discount
profit(x_grid, h_grid[i]) * (1 - delta)
})
# find havest, h that gives the maximum value
out <- sapply(1:gridsize, function(j){
value <- max(V1[j,], na.rm = T) # each col is a diff h, max over these
index <- which.max(V1[j,]) # store index so we can recover h's
c(value, index) # returns both profit value & index of optimal h.
})
# Sets V[t+1] = max_h V[t] at each possible state value, x
V <- out[1,] # The new value-to-go
D[,OptTime-time+1] <- out[2,] # The index positions
}
# Reed derives a const escapement policy saying to fish the pop down to
# the largest population for which you shouldn't harvest:
ReedThreshold <- x_grid[max(which(D[,1] == 1))]
# Format the output
list(D=D, V=V, S=ReedThreshold)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.