simple_spot_rate_extraction.m

HERE IS A GOOD PSO TOOL BOX

% Simplified spot rate curve extraction (coupon stripping).
%
% This is a simplified example from the paper:
%
%   Analyzing the Indonesian Government Bond Market
%   with a New Coupon Stripping Model and Principal Components
%   by K. O. Kortanek and H. Yunianto
%   (see page 15 for the GP problem formulation)
%
% Performs the spot rate curve extraction (coupon stripping) from
% financial data that only include bonds with no coupon payments.
% This results in a GP that can be found on the page 15 of the paper.
%
% Almir Mutapcic 01/15/06

%********************************************************************
% problem data
%********************************************************************
% problem sizes
N = 2;   % number of bonds (length of PVhat)
L = 1;   % number of uniform perturbation intervals

% observed prices based on hypothetical bond data (this is made up data)
PVhat = [98.85; 96.25];
% maturity for each bond (in years)
T = [0.008; 0.015];

% problem constants (all made up)
RH    = 1; RL    = .5;
RHinf = 1; RLinf = .5;
WH    = 1; WL    = .5;

global beta gamma1 gamma2;
beta   = -0.1005;
gamma1 = 0.999*ones(N,1);
gamma2 = 0.999*ones(N,1);

%********************************************************************
% problem variables, objective function, and constraints
%********************************************************************
% GP variables
gpvar x(2*L+2+1)     % note: w_bar = exp(x_2*L+2+1)

% objective function
obj = inv(x(L+2+1)) + x(L+2+1);

% model prices
% denote P_k(T_k) as PV
PV = posynomial;
for k = 1:N
  PV(k,1) = 100*x(1)^ar(T(k))*x(2)^aalpha(T(k));
  PV(k,1) = PV(k)*x(3)^b(T(k)); % here b is the same as aalpha
end

% hard constraint on r
c1 = [ exp(-RH)*x(1) <= 1; exp(RL)*x(1)^(-1) <= 1];

% hard constraint on alpha/-beta
c2 = [ exp(RHinf*beta)*x(2) <= 1; exp(-RLinf*beta)*x(2)^(-1) <= 1];

% "driving" constraints
c_drive = [ exp(-WH)*x(3)*x(4)^-1*x(5)^-1 <= 1;
            exp(WL)*x(3)^-1*x(4)^-1*x(5)  <= 1 ];

% constraint set 
constr = [ c_drive; c1; c2; PV <= PVhat ];

%********************************************************************
% solving the GP problem
%********************************************************************
% extract the spot rate curve from the observables
[oval sol status] = gpsolve(obj, constr);

% evaluate and display computed prices
disp(' ')
fprintf(1,'Estimate of PV(1) is %3.4f versus observed PV_hat(1) = %3.4f.\n',...
        eval(PV(1), sol), PVhat(1));
fprintf(1,'Estimate of PV(2) is %3.4f versus observed PV_hat(2) = %3.4f.\n',...
        eval(PV(2), sol), PVhat(2));