Matlab Calling Functions from Mupad Notebooks in MatlabR2015a
1 次查看(过去 30 天)
显示 更早的评论
After loading "main.m" and the MuPAD notebook, here is the error message I got:
Warning: The system is rank-deficient. Solution is not unique.
> In symengine (line 57)
In sym/privBinaryOp (line 903)
In / (line 297)
In main (line 93)
Error using lsqfcnchk (line 108)
If FUN is a MATLAB object, it must have an feval method.
Error in lsqnsetup (line 37)
funfcn = lsqfcnchk(FUN,caller,lengthVarargin,funValCheck,flags.grad);
Error in lsqnonlin (line 182)
[funfcn,mtxmpy,flags,sizes,funValCheck,xstart,lb,ub,EXITFLAG,Resnorm,FVAL,LAMBDA, ...
Error in main (line 101)
X = lsqnonlin(Call(kappa,theta,delta,vega,rho) - Market,X0,[0,0,0,0,0],[],options);
Actually I do not plan to solve functions l(n),m(n),and c(n) explicitly. They are defined recurrently. How should I fix this bug? Can my idea be fulfilled by Matlab? (maybe Matlab would not do the recurrence equations for me?)
Excel is attached to this question.
main.m
% This Task Consists 2 Steps.
% Step1: Debug "main.m" so that for given values of (K,t,s), "Call" can be expressed as a function of 5 parameters:
% (kappa,theta,delta,vega,rho), namely Call = function(kappa,theta,delta,vega,rho). This expression can be implicit.
% Step2: After reading 100+ sets of numbers (K,t,s) from excel all at once, we get an array of Call's.
% Calibrate the array of Call's to the array of "Market Values" in excel using the Levenberg-Marquardt algorithm,
% a built-in function in Matlab, to get the optimal values of parameters (kappa,theta,delta,vega,rho).
syms kappa theta delta vega rho % 5 parameters
syms tau n
% Case for 2017.10.4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K = xlsread('2017.10.4.xls','B5:B190');
t = xlsread('2017.10.4.xls','C5:C190');
s = 0;
Future = 11.35;
r = 0.0108;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t0 = 0.00001; % Fixed data
g = 1/2; % Fixed data
% Define variables using (kappa,theta,delta,vega,rho) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = 2 * kappa / delta^2 + 1;
betaa = 2 * kappa * theta / delta^2 + 1;
C = 1 / sqrt(2 * pi * vega);
lambda = 1 / (2 * vega);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define function l(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms l(n)
% result = feval(symengine,l,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define nested functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% f1 is a function of tau, which disappears in f2
f1 = C * rho * (1 + t0)^(g * rho - 1) * tau^(-g - 1) * exp(-lambda * tau * (1 + t0)^(-rho)) * (g + lambda * tau * (1 + t0)^(-rho));
% f2 depends on f1; f2 is an integral of tau
f2 = int((1 - exp(-n * kappa * theta * tau)) * f1,tau,0,inf);
% f3 depends on f2 and l(n)
f3 = exp(-t * f2) * sqrt(factorial(n) / gamma(n + alpha)) * l(n);
% A depends on f3; n disappears in A
A = betaa * gamma(alpha - 1) * symsum(f3,n,0,inf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k = K .* A / Future;
% Define function m(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms m(n)
% result = feval(symengine,m,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define function c(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
syms c(n)
% result = feval(symengine,c,n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define function B
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B = exp((s - t) * f2) * c(n) * l(n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Express Call implicitly using (kappa,theta,delta,vega,rho)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = [kappa,theta,delta,vega,rho];
Call = symfun(exp(-r * t) * Future / A * symsum(B,n,0,inf), X);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform the Levenberg-Marquardt algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X0 = [0,0,0,0,0];
Market = xlsread('2017.10.4.xls','F5:F190');
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
X = lsqnonlin(Call(kappa,theta,delta,vega,rho) - Market,X0,[0,0,0,0,0],[],options);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Here is the code in a single MuPAD notebook(hit "Enter" after typing each equation):
solve(rec(l(n) = (alpha - 2 + 2 * n - betaa) / sqrt(n * (alpha - 1 + n)) * l(n - 1) - sqrt((alpha + n - 2) * (n - 1) / (alpha - 1 + n) / n) * l(n - 2),l(n),{l(0) = 1 / sqrt(gamma(alpha)), l(1) = (alpha - betaa) / sqrt(gamma(alpha + 1))}))
solve(rec(m(n) = (alpha + 2 * n - 1 - betaa / k) / sqrt(n * (alpha + n)) * m(n - 1) - sqrt((alpha + n - 1) * (n - 1) / (alpha + n) / n) * m(n - 2),m(n),{m(0) = 1 / sqrt(gamma(alpha + 1)), m(1) = (1 + alpha - betaa / k) / sqrt(gamma(alpha + 2))}))
solve(rec(c(n) = c(n-1) * sqrt(n / (alpha + n - 1)) + betaa^alpha / k^(alpha - 1) * m(n - 2) * exp(-betaa / k) / sqrt(n * (n - 1) * (n + alpha - 1)),c(n),{c(0) = sqrt(gamma(alpha)) * (betaa / (alpha - 1) - k) * gammainc(betaa / k,alpha - 1) + (betaa^(alpha - 1) / k^(alpha - 2)) * exp(-betaa / k) / sqrt(gamma(alpha)),c(1) = sqrt(gamma(alpha + 1)) * betaa / alpha / (alpha - 1) * (1 - igamma(alpha - 1,betaa / k) / gamma(alpha - 1))}))
0 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!