Returning multiple values in fminbnd

3 views (last 30 days)
Prerna Mishra
Prerna Mishra on 25 Feb 2023
Commented: Torsten on 25 Feb 2023
I am trying to solve a dynamic optimization problem in which I solve for optimxal capital and debt. The problem is structured such that for each period I have to find the value of debt, variable d, that maximizes the firm value, given the value of k and k0. And then in the next step I have to find the value of k that maximizes the firm value over time, given k0. My main looks as follows:
disc = 0.96; %discount rate for the expected next period value function
alpha = 0.45; %capital elasticity, Cobb-Douglas *
beta = 0.34; %labor elasticity, Cobb-Douglas *
R = 719.0615; %Market Size
f = 4.02; %fixed cost of production
sigma = 5.23; %elasticity of substitution
W = 0.3; %wage *
theta = 0.7; %cost of bankruptcy *
chi = 0.8; %level of debt at which bankruptcy happens *
c = 0.2; %parameters of capital adjustment *
kappa = 0.3; %managerial rent *
mu_b = 0.1; %mean of the firm level shock *
sigma_b = 0.12; %standard deviation of the firm level shock
omega = (beta*(sigma-1))/(beta*(sigma-1) - sigma);%an itermediary parameter
na = 3;%TFP shock grid size
nk = 11; %capital grid size
%setting up the capital grid
kmin = 1;
kmax = 10000000000;
kvec = linspace(log(kmin),log(kmax),nk)';
kvec = exp(kvec);
%setting up the TFP grid
multiple = 3;
mu_a = 0;
rho_a = 0.7;
sigma_a = 0.05;
[avec,pi]=tauchen(0.0,sigma_a,rho_a,multiple,na);
Unrecognized function or variable 'tauchen'.
avec = exp(avec)';
tol = 1e-5;
diff = 1;
its = 0;
maxits = 500;
v0 = zeros(nk,na);
while diff>tol && its <maxits
for j = 1:na
for i = 1:nk
P0 = 0.001;
a0 = avec(j);
k0 = kvec(i);
k1 = fminbnd(@(k)ValueFunctionCapital(k,v0,k0,a0,j,P0),kmin,kmax);
[val,d] = -ValueFunctionCapital(k1,v0,k0,a0,j,P0);
P = fixedPointP(P0,k1);
v1(i,j) = val;
k11(i,j) = k1;
d11(i,j) = d;
p11(i,j) = P;
end
end
diff = norm(v0-v1)
v0 = v1;
its = its+1
end
My ValueFuntionFunctionCapital looks as follows:
function [val,d] = ValueFunctionCapital(k,v0,k0,a0,j,P)
global disc alpha beta R f sigma W theta chi c kappa mu_b sigma_b omega pi mu_a sigma_a kvec avec maxits
g1 = interp1(kvec,v0,k,'linear'); %smooth out previous value function
%some constants
prof_brac = (sigma/(beta*(sigma-1)))^omega - (sigma/(beta*(sigma-1)))...
^(sigma/(beta*(sigma-1)-sigma));
quad_cost = c*(k0-k)^2;
function val_inner = ValueFunctionDebtInner(d,k0,k,a0,P)
g_inner = @(b) lognpdf(b, mu_b, sigma_b); % define PDF of lognormal distribution
fun_inner = @(b) (1./b).^omega .* g_inner(b);
Q11_inner = integral(fun_inner, 0, Inf);
pdf_inner = @(x) lognpdf(x, mu_b, sigma_b);
V_term1_1_inner = (1-kappa)*(R*P^(sigma-1))^(-1/(beta*(sigma-1)))*(a0*k^alpha)^...
(-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*prof_brac*Q11_inner;
V_term1_2_inner = (1-kappa)*(k0-k-f*k-quad_cost);
V_term1_inner = V_term1_1_inner + V_term1_2_inner;
b_th_inner = bankruptcy(a0,k0,k,d,P,prof_brac,W,chi,kappa,f,sigma,beta,c,R);
V_term2_cdf_inner = logncdf(b_th_inner, mu_b, sigma_b); % compute CDF at x ;
V_term2_inner = kappa*d*(1-V_term2_cdf_inner);
Q31_inner = integral(fun_inner, 0, b_th_inner);
% define PDF function using lognpdf
Q32_inner = integral(pdf_inner, 0, b_th_inner);
V_term3_1_inner = (1-kappa)*theta*(R*P^(sigma-1))^(-1/(beta*(sigma-1)-sigma))*...
(a0*k^alpha)^ (-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*...
prof_brac*Q31_inner;
V_term3_2_inner = (1-kappa)*theta*(k0-k-f*k-quad_cost)*Q32_inner;
V_term3_inner = V_term3_1_inner + V_term3_2_inner;
V_final_inner = V_term1_inner + V_term2_inner - V_term3_inner ;
val_inner = V_final_inner;
val_inner = -val_inner;
end
dmin = 1;
dmax = k;
d = fminbnd(@(d)ValueFunctionDebtInner(d,k0,k,a0,P),dmin,dmax);
g = @(b) lognpdf(b, mu_b, sigma_b); % define PDF of lognormal distribution
fun = @(b) (1./b).^omega .* g(b);
Q11 = integral(fun, 0, Inf);
pdf = @(x) lognpdf(x, mu_b, sigma_b);
V_term1_1 = (1-kappa)*(R*P^(sigma-1))^(-1/(beta*(sigma-1)))*(a0*k^alpha)^...
(-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*prof_brac*Q11;
V_term1_2 = (1-kappa)*(k0-k-f*k-quad_cost);
V_term1 = V_term1_1 + V_term1_2;
b_th = bankruptcy(a0,k0,k,d,P,prof_brac,W,chi,kappa,f,sigma,beta,c,R);
V_term2_cdf = logncdf(b_th, mu_b, sigma_b); % compute CDF at x ;
V_term2 = kappa*d*(1-V_term2_cdf);
Q31 = integral(fun, 0, b_th);
% define PDF function using lognpdf
Q32 = integral(pdf, 0, b_th);
V_term3_1 = (1-kappa)*theta*(R*P^(sigma-1))^(-1/(beta*(sigma-1)-sigma))*...
(a0*k^alpha)^ (-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*...
prof_brac*Q31;
V_term3_2 = (1-kappa)*theta*(k0-k-f*k-quad_cost)*Q32;
V_term3 = V_term3_1 + V_term3_2;
%V_term3 = -V_term3;
V_final = V_term1 + V_term2 - V_term3 ;
val = V_final + disc * (g1*pi(j,:)');
val = -val;
end
I understand that there is a lot of repetition of code here. But for now what I want is that I can store k, the value of capital that maximizes the value and the corresponding d. I am having troubles storing that since with fminbnd I am not able to retrieve both the values when I call the function. Is there anyway I can get both k and corresponding d in main?
  3 Comments
Prerna Mishra
Prerna Mishra on 25 Feb 2023
k1 = fminbnd(@(k)ValueFunctionCapital(k,v0,k0,a0,j,P0),kmin,kmax);
In the above line I want both the optimized k and d.
Torsten
Torsten on 25 Feb 2023
You can get d only by calling "ValueFunctionCapital" again with the optimal k1 as you do in the next line of your code. Doesn't that work ? I cannot run your code since at least "tauchen" is missing.

Sign in to comment.

Answers (0)

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!