Algorithms for ‘Selling to the Newsvendor: An analysis of Price-only Contracts’
I developed the algorithms during my PhD study to deepen my understanding of both MATLAB and Newsvendor model. I found it was quite useful. Here I also put the algorithms for your reference. It would be appreciated if you can send your comments when spotting any errors.
Analysis 1
clc;
syms y a b c r
% specify variables
unifcdf = @unifCDF;
% call the Uniform Cumulative Distribution Function (CDF)
optiFun_y = (1 - unifcdf(y, a, b)) * (1 - y * diff(unifcdf(y, a, b),y) / (1 - unifcdf(y, a, b))) - c/r
% the equation used to calcualte optimal order quantity
%%
clc;
syms y
A = [100, 100, 180];
B = [500, 600, 500];
% parameters regarding different market sizes
M = zeros(3, 4);
% create an empty matrix to store parameters and results
for i = 1 : 3
x = [4, 10, A(i), B(i)];
c = x(1);
r = x(2);
a = x(3);
b = x(4);
opti_unif_y = solve(((a - y)/(a - b) - 1)*(y/(((a - y)/(a - b) - 1)*(a - b)) - 1) - c/r);
% calculate the optimal order quantity
% the equation in the solve function is produced by 'optiFun_y' in the
% last section
opti_unif_w = double(r * (1 - unifcdf(opti_unif_y, a, b)));
% calculate the optimal wholesale price
M(i, 1 : 4) = [a, b, opti_unif_y, opti_unif_w];
% grap the data to the empty matrix
end
M
% show the matrix with stored values
%%
%################### functions for numerical analysis 1 ##################%
%=================== set up the function of the uniform CDF ==============%
function uniCdf = unifCDF(x, a, b)
uniCdf = (x -a) / (b - a);
end
endAnalysis 2
clc;
syms y k theta
% specify variables
assume(k > 0 & theta >0)
% set assumption on symbolic object
c = 4;
r = 10;
% assign certain parameters values
powcdf = @powCDF;
% call the Power Cumulative Distribution Function (CDF)
%=== calculate the optimial order in the decentralised supply chain (DSC) ====%
optiPowFun_y_dsc = (1 - powcdf(y, theta, k)) * (1 - y * ...
diff(powcdf(y, theta, k),y) / (1 - powcdf(y, theta, k))) - c/r
% this function is for calculating the optimal quantity in the DSC
optiQuan_ds = solve(optiPowFun_y_dsc, y)
% the optimal quantity for DSC
%===== calculate the optimial order in the integrated supply chain (ISC) ====%
optiPowFun_y_isc = powcdf(y, theta, k) - 1 + c/r;
% the function for calculating the optimal quantity in ISC
optiQuan_isc = solve(optiPowFun_y_isc, y)
% the optimal quantity for ISC
%=============== define coefficient of variation (CV) ===================%
syms y
% specify the variable
expect_x = int(y * diff(powcdf(y, theta, k), y), 0, 1);
% calculate expected value
sd_x = sqrt(int(y^2 * diff(powcdf(y, theta, k), y), 0, 1) - expect_x^2);
% calculate standard deviation
cvFun = sd_x / expect_x
% calculate CV
%%
clc;
syms y
powcdf = @powCDF;
% call the Power CDF
K = linspace(0.2, 8, 30);
% generate linearly spaced vector for k
for i = 1 : 30
c = 4;
r = 10;
theta = 1;
k = K(i)
% assign certain parameters values
cv(i) = -(1 - (1/(k + 1) - 1)^2 - 2/(k + 2))^(1/2)/(1/(k + 1) - 1);
% coefficient of variation (CV)
optimal_y = double(theta*(3/(5*(k + 1)))^(1/k));
% the optimal order quantity
prop_dm = double(((1 - powcdf(optimal_y, theta, k)) - c/r) * optimal_y);
% the numerator of the ratio (the manufacturer's profit)
prop_dsc = double(optimal_y * (1 - powcdf(optimal_y, theta, k)) +...
int(y * diff(powcdf(y, theta, k), y), 0, optimal_y) - c/r * optimal_y);
% the denominator of the ratio (the DSC's profit)
prop_dm_dsc(i) = double(prop_dm / prop_dsc)
% the ratio of the manufacturer's to DSC's profit
optiQuan_is(i) = (3/5)^(1/k)*theta;
% the optimal quanity of ISC
prop_isc = - c/r * optiQuan_is(i) + optiQuan_is(i) -...
int(powcdf(y, theta, k), 0, optiQuan_is(i));
% denominator of the ratio (ISC's profit)
prop_dsc_isc(i) = double(prop_dsc/prop_isc)
% the ratio of the decentralised to ISC profit
end
%%
plot(K, prop_dm_dsc * 100, 'k-.')
% plot the ratio of the manufacturer's to DSC's profit
ytickformat('percentage')
% specify y-axis tick label format
hold on;
% retains the current plot and certain axes properties
% so that subsequent graphing commands add to the existing graph.
[hAx, hLine1, hLine2] = plotyy(K, prop_dsc_isc * 100, K, cv)
% plot the ratio of the decentralized to integrated supply chain profit
% plot the relationship between K and cv
title('Supply-Chain Performance Under Power Function Demand')
xlabel('k')
ylabel(hAx(2), 'Coefficient of Variation') % right y-axis
% add the title and name the x and y axis
set(hLine1, 'Color', 'k', 'LineStyle', '-');
set(hLine2, 'Color', 'k', 'LineStyle', '-', 'Marker','X');
% set the color, style and marker of the line
set(hAx, 'XTick', 0: 2 : 8, {'YTick'}, {0: 10 :100; 0 : 0.25 : 2}, {'ycolor'},{'k';'k'})
% customising the tick values
legend('\Pi_M/\Pi_D', '\Pi_D/\Pi_I','CV','Location','east')
% add legend
box off;
% removes the box outline around the current axes
hold off;
%%
%#################### functions for numercial analysis 2 #################%
%=================== CDF of the power distribution =======================%
function pwercdf = powCDF(y, theta, k)
pwercdf = (y / theta)^k;
endAnalysis 3
%##################### numerical analysis 3 ##############################%
clc;
syms y k
% specify variables
M_num1 = zeros(4, 8);
% create an empty matrix to store parameters and results
for j = 1 : 4
K = [0.20, 0.25, 0.50, 1.00, 2.00, 4.00, 9.00, 16.00];
C = [2, 4, 6, 8];
c = C(1);
r = 10;
theta = 1.2;
k = K(j)
% assign certain parameters values
gamapdf = @gamaPDF;
% call the Gamma Cumulative Distribution Function (CDF)
expect_gam_x = int(y * gamapdf(y, theta, k), y, 0, +inf);
sd_gam_x = sqrt(int(y^2 * gamapdf(y, theta, k), y, 0, +inf) - expect_gam_x^2);
cv_gam = double(sd_gam_x / expect_gam_x);
% define coefficient of variation (CV)
%========= calculate the optimial order in the decentralized supply chain (DSC) ====%
optiGamFun_y_dsc = (1 - int(gamapdf(y, theta, k),0 , y)) * (1 - y *...
gamapdf(y, theta, k) / (1 - int(gamapdf(y, theta, k),0 , y))) - c/r;
% the function for calculating the optimal quantity of DSC
optiQuan_ds_gama = double(vpasolve(optiGamFun_y_dsc, y));
% the optimal quanity for DSC
opti_w_gama = r * (1 - int(gamapdf(y, theta, k),0, optiQuan_ds_gama));
% the optimal wholesale price
ratio_w_r = double(1 - int(gamapdf(y, theta, k),0, optiQuan_ds_gama));
% the tatio of wholesale price to retail price
%========= calculate optimial order in the integrated supply chain =======%
optiGamaFun_y_isc = int(gamapdf(y, theta, k),0, y) + c/r - 1;
% the function used to calculate the optimal quantity of integrated
% supply chain (ISC)
optiQuan_isc_gama = double(vpasolve(optiGamaFun_y_isc, y));
% the optimal quanity for ISC
ratio_ydsc_yisc = optiQuan_ds_gama / optiQuan_isc_gama;
% the ratio of DSC to ISC order quantity
prop_dsc_profit = (1 - int(gamapdf(y, theta, k),0, optiQuan_ds_gama)) *...
optiQuan_ds_gama + int(y * gamapdf(y, theta, k), 0, optiQuan_ds_gama)...
- c/r * optiQuan_ds_gama;
% the numerator of the ratio (the DSC's profit)
prop_isc_profit = - c/r * optiQuan_isc_gama + optiQuan_isc_gama -...
int(int(gamapdf(y, theta, k),0, y), 0, optiQuan_isc_gama);
% the denominator of the ratio (the ISC's profit)
ratio_dsc_isc_profit = double(prop_dsc_profit / prop_isc_profit);
% the ratio of DSC to the ISC profit
prop_m_profit = (1 - int(gamapdf(y, theta, k),0, optiQuan_ds_gama) - c/r) * optiQuan_ds_gama;
ratio_m_dsc_profit = double(prop_m_profit / prop_dsc_profit);
% the ratio of manufactuer to DSC profit
ratio_m_isc_profit = double(prop_m_profit / prop_isc_profit);
% the ratio of manufactuer to ISC profit
M_num1(j, 1:8) = [c/r, k, cv_gam, ratio_ydsc_yisc * 100,...
ratio_w_r, ratio_dsc_isc_profit*100, ratio_m_dsc_profit*100,...
ratio_m_isc_profit*100]
% grap the data to the empty matrix
end
M_num1
% show the matrix with stored values
%%
%#################### functions for numercial analysis 3 #################%
%============================== define Gamma CDF =========================%
function gamapdf = gamaPDF(y, theta, k)
gamapdf = theta^(-k) * y^(k -1) * exp(-y/theta)/gamma(k);
endAnalysis 4
clc;
syms y k
% specify variables
assume(k > 0)
% set assumption on symbolic object
c = 4;
r = 10;
theta = 1;
% assign certain parameters values
powcdf = @powCDF;
% call the Power Cumulative Distribution Function (CDF)
%=== calculate the optimial order in the decentralised supply chain (DSC) ====%
optiPowFun_y_dsc = (1 - powcdf(y, theta, k)) * (1 - y * ...
diff(powcdf(y, theta, k),y) / (1 - powcdf(y, theta, k))) - c/r;
% this function is for calculating the optimal quantity in the DSC
optiQuan_ds = solve(optiPowFun_y_dsc, y)
% the optimal quantity for DSC
%===== calculate the optimial order in the integrated supply chain (ISC) ====%
optiPowFun_y_isc = powcdf(y, theta, k) - 1 + c/r;
% the function for calculating the optimal quantity in ISC
optiQuan_isc = solve(optiPowFun_y_isc, y)
% the optimal quantity for ISC
%=============== define coefficient of variation (CV) ===================%
syms y
% specify the variable
expect_x = int(y * diff(powcdf(y, theta, k), y), 0, 1);
% calculate expected value
sd_x = sqrt(int(y^2 * diff(powcdf(y, theta, k), y), 0, 1) - expect_x^2);
% calculate standard deviation
cvFun = sd_x / expect_x
% calculate CV
%%
clc;
syms y
assume(y>0)
powcdf = @powCDF;
% call the Power CDF
M_num3 = zeros(32, 7)
% create an empty matrix to store parameters and results
for j = 1 : 4
BETA = [0.3, 0.4, 0.5, 0.6];
for i = 1 : 8
K = [0.1, 0.25, 0.414, 0.75, 1.00, 1.5, 2.00, 3.00];
c = 4;
r = 10;
theta = 1;
k = K(i)
% assign certain parameters values
cv = -(1 - (1/(k + 1) - 1)^2 - 2/(k + 2))^(1/2)/(1/(k + 1) - 1);
% coefficient of variation (CV)
optimal_y = double(solve(- ((k*y*y^(k - 1))/(y^k - 1) + 1)*(y^k - 1) - 2/5, y));
% the optimal order quantity
prop_dm =((1 - powcdf(optimal_y, theta, k)) - c/r) * optimal_y;
% the numerator of the ratio (the manufacturer's profit)
prop_dsc = optimal_y * (1 - powcdf(optimal_y, theta, k)) +...
int(y * diff(powcdf(y, theta, k), y), 0, optimal_y) - c/r * optimal_y;
% the denominator of the ratio (the DSC's profit)
prop_dm_dsc = prop_dm / prop_dsc;
% the ratio of the manufacturer's to DSC's profit
optiQuan_is = (3/5)^(1/k);
% the optimal quanity of IDS
prop_isc = - c/r * optiQuan_is + optiQuan_is - int(powcdf(y, theta, k), 0, optiQuan_is);
% denominator of the ratio (ISC's profit)
prop_dsc_isc = prop_dsc/prop_isc;
% the ratio of the DSC to ISC profit
%========== optimal strategy considering the retailer's power=====%
optiQuanBe = @optiQuanBeta;
beta = BETA(j)
% assign certain parameters values
optiQuan_beta = optiQuanBe(c, r, beta, k);
% the optimal order quanity considering the retailer's power
optiQuan_beta_SC = optiQuan_beta / optimal_y;
% ratio of beta SC to integrated SC quantity
prop_dsc_beta = optiQuan_beta * (1 - powcdf(optiQuan_beta, theta, k)) +...
int(y * diff(powcdf(y, theta, k), y), 0, optiQuan_beta) - c/r * optiQuan_beta;
% denominator of the ratio (beta SC's profit)
ratio_betaSC_isc = double(prop_dsc_beta / prop_isc);
% ratio of beta SC to ISC profit
ratio_betaSC_dsc = double(prop_dsc_beta/prop_dsc);
% ratio of beta SC to DSC profit
prop_beta_m = ((1 - powcdf(optiQuan_beta, theta, k)) - c/r) * optiQuan_beta;
% the numerator of the ratio (the manufacturer's profit)
ratio_betaM_DscM = prop_beta_m / prop_dm;
% ratio of the manufacturer's profit in beta SC to the
% manufacturer's profit in DSC.
M_num3(8 * (j-1)+ i, 1:7) = [beta, k, cv, optiQuan_beta_SC*100,...
ratio_betaSC_isc*100, ratio_betaSC_dsc*100, ratio_betaM_DscM*100];
% grap the data to the empty matrix
end
end
M_num3
% show the matrix with stored values
%%
%#################### functions for numercial analysis 4 #################%
%================== power CDF ========================%
function pwercdf = powCDF(y, theta, k)
pwercdf = (y / theta)^k;
end
%============ the optimal quantity considering the retailer' power ========%
function optiQuan = optiQuanBeta(c, r, beta, k)
optiQuan = ((r - c) * beta * (k + 1) / (r * (k + beta)) )^(1/k);
end