56 lines
1.8 KiB
R
56 lines
1.8 KiB
R
|
# Transition dynamics
|
|||
|
|
# install.packages("pracma")
|
|||
|
|
library(pracma)
|
|||
|
|
|
|||
|
|
## For forward shooting
|
|||
|
|
a1 <- function(t){
|
|||
|
|
term_a1 = λ_1 *
|
|||
|
|
((1- λ_1)/( λ_1 * w_3 * exp(γ_1 * t)))^((1-λ_1) * (1-η)) /
|
|||
|
|
(k1_forward_shooting * exp((ρ-r_1) * t))
|
|||
|
|
result_a1 = term_a1^(1/η)
|
|||
|
|
return(result_a1)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
a3 <- function(t){
|
|||
|
|
first_term_a3 = w_3 * exp(γ_1 * t) * δ_3 * λ_3 *
|
|||
|
|
β_3 * exp(γ_3 * t) / (r_1 - γ_1 - r_3)
|
|||
|
|
second_term_a3 = (1-λ_3)/(λ_3 * w_3 * exp(γ_1 * t))
|
|||
|
|
subresult_a3 = first_term_a3 * (second_term_a3^(δ_3*(1-λ_3)))
|
|||
|
|
result_a3 = subresult_a3^(1/(1-δ_3))
|
|||
|
|
return(result_a3)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
s1 <- function(t){
|
|||
|
|
result_s1 = ((1-λ_1)/ λ_1) * (a1_t / (x3_t * w_3 * exp(γ_1 * t)))
|
|||
|
|
return(result_s1)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
s3 <- function(t){
|
|||
|
|
result_s3 = ((1-λ_3)/ λ_3) * (a3_t / (x3_t * w_3 * exp(γ_1 * t)))
|
|||
|
|
return(result_s3)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
dx1 <- function(t){
|
|||
|
|
result_x1 = - a1_t - a3_t + x3_t * w_3 * exp(γ_1 * t) * (1-s1_t - s3_t)
|
|||
|
|
return(result_x1)
|
|||
|
|
}
|
|||
|
|
dx3 <- function(t){
|
|||
|
|
result_x3 = β_3 * exp(γ_3 * t) * (a3(t)^λ_3 * (s3_t*x3_t)^(1-λ_3))^δ_3
|
|||
|
|
return(result_x3)
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
dx3_reverse_shooting <- function(t, guessmin, guessmax, x3_target, t_previous){
|
|||
|
|
## Inverting dx3 turns out to be a little bit involved.
|
|||
|
|
check_hypothesis <- function(x3_hypothesis){
|
|||
|
|
dx3_hypothesis = β_3 * exp(γ_3 * t_previous) * (a3(t_previous)^λ_3 * (s3_t*x3_hypothesis)^(1-λ_3))^δ_3
|
|||
|
|
difference = x3_target - (x3_hypothesis + dx3_hypothesis*stepsize)
|
|||
|
|
return(difference)
|
|||
|
|
}
|
|||
|
|
x3_best_guess <- bisect(check_hypothesis, guessmin, guessmax)$root
|
|||
|
|
#uniroot(check_hypothesis, lower=guessmin, upper=guessmax, tol = 10^(-22))[[1]]
|
|||
|
|
# x3_best_guess <- uniroot(check_hypothesis, lower=(x3_target-2), upper=(x3_target), tol = 10^(-22), extendInt = "down")[[1]]
|
|||
|
|
## x3_best_guess <- uniroot(check_hypothesis, c(0, x3_target), tol = 10^(-22), a = 0, extendInt = "upX")[[1]]
|
|||
|
|
return(x3_best_guess)
|
|||
|
|
}
|
|||
|
|
|