Subscript indices must either be real positive integers or logicals.

function [r,fun_val]=gradient_method_backtrackingMbb(r0,s,alapha,beeta,epsilon)

% r=[0.2 0.3 0.2 0.4 0.2 0.5 0.2 0.03 0.3 0.3 0.3 0.2];%% Risk-Free Interest Rate

sigma=0.5; %% Target Volatility from the market

X=100; %% Example Strike Price

T=1; %% Years until Expiry

%%

%% Grid Parameters %%

M=11; %% Asset Grid points

N=11; %% Time Grid points

xzero=0; %% Specify Minimum Share Price

Smax=200; %% Specify Maximum Share Price

%%

%% Grid Setup and Boundary Conditions %%

j=0:M; %% Set up j vector

dt=T/N; %% Construct time step

ds=Smax/M; %% Construct price step

ug=zeros(M+1,N+1); %% Initialise Grid i.e. matrix

xgrid=0:ds:Smax; %% Positive, equally spaced price steps. %% i.e. elements

Tgrid=T:-dt:0; %% Negative(Backward), equally spaced ti %% time steps

%% Boundary Conditions

ug(1,:)=X*exp(-r0.*Tgrid); %% Boundary Condition : Price=0

ug(:,end)=max(X-xgrid,0); %% Boundary Condition : Payoff of Put

ug(end,:)=0; %% Boundary Condition : Price tending to %% to "infinity"

%%

alpha=-(1/2)*dt*(sigma^2*j.^2); % function alpha

beta=1+(sigma^2*j.^2 +r0.*j+r0.*j)*dt; % function beta

gamma=-(1/2)*dt*(r0.*j+sigma^2*j.^2); % function gamma

%% Construction of Soln Matrix %%

A=diag(alpha(3:M),-1)+diag(beta(2:M))+diag(gamma(2:M-1),1);

% Here we create Matrix (33) with betas on the leading diagonal, alphas

% offset down(-1) and gammas offset up (+1)

Ainv=inv(A); % Create inverse of A to test stability

Anormi=norm(Ainv,inf); % Stability Test

C=zeros(size(A,2),1); %Create matrix of 0's w/1 column and number of rows= =no. of columns in A

for i=N:-1:1 %For loop solves our M-1 eqns for every time grid point N=2310

C(1)=alpha(2)*ug(1,i); % first element of C

C(end) = gamma(end)*ug(end,i); % Will always be zero as previously stated hence irrelevant

ug(2:M,i)=A\(ug(2:M,i+1)-C); %Inverted matrix soln for P_i

end

%Approximate Theoretical Data

% xx=138.50; %% Initial Share Price

% r=[0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01]; %% Risk-Free Interest Rate

sigmak=[0.100000000000000,0.156346511368286,0.208128163491120,0.251149914870852,0.281926399070904,0.297964288376187,0.297964288376187,0.281926399070904,0.251149914870852,0.208128163491120,0.156346511368286,0.100000000000000]; %%01 Volatility

X=100; %% Example Strike Price

T=1; %% Years until Expiry

%%

%% Grid Parameters %%

M=11; %% Asset Grid points

N=11; %% Time Grid points

% xxzero=0; %% Specify Minimum Share Price

Smax=200; %% Specify Maximum Share Price

%%

%% Grid Setup and Boundary Conditions %%

j=0:M; %% Set up j vector

dt=T/N; %% Construct time step

ds=Smax/M; %% Construct price step

u=zeros(M+1,N+1); %% Initialise Grid i.e. matrix

xgrid=0:ds:Smax; %% Positive, equally spaced price steps. %% i.e. elements

Tgrid=T:-dt:0; %% Negative(Backward), equally spaced ti %% time steps

%% Boundary Conditions

u(1,:)=X*exp(-r0.*Tgrid); %% Boundary Condition : Price=0

u(:,end)=max(X-xgrid,0); %% Boundary Condition : Payoff of Put

u(end,:)=0; %% Boundary Condition : Price tending to %% to "infinity"

alpha1=-(1/2).*dt.*(sigmak.^2.*j.^2); % function alpha

beta1=1+(sigmak.^2.*j.^2 +r0.*j+r0.*j)*dt; % function beta

gamma1=-(1/2).*dt.*(r0.*j+sigmak.^2.*j.^2); % function gamma

%% Construction of Soln Matrix %%

B=diag(alpha1(3:M),-1)+diag(beta1(2:M))+diag(gamma1(2:M-1),1);

% Here we create Matrix (33) with betas on the leading diagonal, alphas

% offset down(-1) and gammas offset up (+1)

Binv=inv(B); % Create inverse of A to test stability

Bnormi=norm(Binv,inf); % Stability Test

C=zeros(size(B,2),1); %Create matrix of 0's w/1 column and number of rows= =no. of columns in C

for i=N:-1:1 %For loop solves our M-1 eqns for every time grid point N

C(1)=alpha1(2)*u(1,i); % first element of C

C(end) = gamma1(end)*u(end,i); % Will always be zero as previously previously stated hence irrelevant

u(2:M,i)=B\(u(2:M,i+1)-C); %Inverted matrix soln for P_i

end

%% Theoretical Volatility V Data from equation V

% Observed Market Data

% S=50; %% Initial Share Price

% r=0.01; %% Risk-Free Interest Rate

% sigma=0.3; %% Volatility

X=100; %% Example Strike Price

T=1; %% Years until Expiry

%%

%% Grid Parameters %%

M=11; %% Asset Grid points

N=11; %% Time Grid points

Szero=0; %% Specify Minimum Share Price

Smax=200; %% Specify Maximum Share Price

%%

%% Grid Setup and Boundary Conditions %%

j=0:M; %% Set up j vector

dt=T/N; %% Construct time step

ds=Smax/M; %% Construct price step

V=zeros(M+1,N+1); %% Initialise Grid i.e. matrix

xgrid=0:ds:Smax; %% Positive, equally spaced price steps. %% i.e. elements

Tgrid=T:-dt:0; %% Negative(Backward), equally spaced ti %% time steps

%% Boundary Conditions

V(1,:)=X*exp(-r0.*j.*Tgrid); %% Boundary Condition : Price=0

V(:,end)=max(X-xgrid,0); %% Boundary Condition : Payoff of Put

V(end,:)=0; %% Boundary Condition : Price tending to %% to "infinity"

alpha2=(1/2)*dt.*(sigmak.^2.*j.^2); % function alpha

beta2=1-(sigmak.^2.*j.^2 +r0.*j+r0.*j)*dt; % function beta

gamma2=(1/2).*dt.*(r0.*j+sigmak.^2.*j.^2); % function gamma

%% Construction of Soln Matrix %%

C=diag(alpha2(3:M),-1)+diag(beta2(2:M))+diag(gamma2(2:M-1),1);

% Here we create Matrix with betas on the leading diagonal, alphas

% offset down(-1) and gammas offset up (+1)

Cinv=inv(C); % Create inverse of A to test stability

Cnormi=norm(Cinv,inf); % Stability Test

Cc=zeros(size(C,2),1); %Create matrix of 0's w/1 column and number of rows= =no. of columns in C

for i=N:-1:1 %For loop solves our M-1 eqns for every time grid point

Cc(1)=alpha2(2)*V(1,i); % first element of C

Cc(end) = gamma2(end)*V(end,i); % Will always be zero as previously previously stated hence irrelevant

V(2:M,i)=C\(V(2:M,i+1)-Cc); %Inverted matrix soln for V_i

end

% Objective function F(x) to minimize in order to solve F(r)=0

F=sum(sum(ug-u).^2);

% Gradient of F (partial derivatives)

dF=(-2*sum(sum(ug-u)*V));

% Gradient method with backtracking stepsize rule

%

% INPUT

%=======================================

% f ......... objective function

% g ......... gradient of the objective function

% x0......... initial point

% s ......... initial choice of stepsize

% alpha ..... tolerance parameter for the stepsize selection

% beta ...... the constant in which the stepsize is multipliedat each backtracking step (0<beta<1)

% epsilon ... tolerance parameter for stopping rule

% OUTPUT

%=======================================

% x ......... optimal solution (up to a tolerance)

% of min f(x)

% fun_val ... optimal function value

r=r0;

grad=dF;

fun_val=F;

iter=0;

while (norm(grad)>epsilon)

iter=iter+1;

t=s;

while (fun_val-F(r-t*grad)<alapha*t*norm(grad)^2)

t=beeta*t;

end

r=r-t*grad;

fun_val=F;

grad=dF;

fprintf('iter_number = %3d norm_grad = %2.6f fun_val = %2.6f \n',...

iter,norm(grad),fun_val);

end