Would anyone have a code (pref. Matlab or R) for any type of estimation (QML, GMM) not using option prices of a stochastic volatility model driven by a CIR process described below?

\begin{equation} dS_t = \mu dt + \sqrt{v_t} dW_t^1 \end{equation} \begin{equation} dv_t = \kappa (\theta - v_t) dt + \sigma \sqrt{v_t} dW^2_t \end{equation} such that $dW_{t}^{1}\,dW_{t}^{2}=\rho dt$

I just need to cross-check my results.


  • $\begingroup$ That's the Heston model. You will find plenty of information on this site and also implementations on Google. Plus, in R it's already implemented (packages fOptions and NMOF). Remember that the Heston model models also the correlation between the two Wiener processes. $\endgroup$
    – Arrigo
    Dec 23, 2014 at 10:52
  • $\begingroup$ @Arrigo: Thanks, there's been a typo in the price process specification. Unlike for Heston, I need to estimate without using the option prices. $\endgroup$ Dec 23, 2014 at 11:16
  • $\begingroup$ Here it is! Heston calibration in Matlab $\endgroup$ Dec 23, 2014 at 12:52
  • $\begingroup$ Answers should contain enough information to actually answer the question. A simple link to another site should just be a comment. $\endgroup$ Dec 23, 2014 at 14:08
  • 1
    $\begingroup$ @GabrielePompa thanks. but this is not a heston model (doesn't assume dependency between the wieners) + for calibration of the heston option prices are needed, which i do not have. $\endgroup$ Dec 23, 2014 at 16:37

1 Answer 1

  1. Purpose:
    • CIR process maximum likelihood estimation.
  2. Input
    • Model.Data = Time series of interest rates observations.
    • Model.TimeStep = Delta t; recommended: 1/250 for daily data and 1/12 for monthly data.
    • Model.Disp = 'y' | 'n', (default: 'n').
    • Model.Method = 'ncx2pdf' | 'besseli' (default: 'besseli').
    • Model.MatlabDisp = 'off' | 'iter' | 'notify' |' final' (default: 'off').
  3. Output:
    • Results.Params = Estimated parameters (kappa, theta, sigma).
    • Results.Fval = Objective function value.

function Results = CIRestimation(Model)
Nobs = length(Model.Data);
r = Model.Data(1:end-1);
dr = diff(Model.Data);           
dr = dr./r.^0.5;
regressors = [Model.TimeStep./r.^0.5, Model.TimeStep*r.^0.5];
drift = regressors\dr;
res = regressors*drift - dr;
alpha = -drift(2);
theta = -drift(1)/drift(2);
sigma = sqrt(var(res, 1)/Model.TimeStep);
InitialParams = [kappa theta sigma];
if ~isfield(Model, 'Disp'), Model.Disp = 'y'; 
if strcmp(Model.Disp, 'y')
    fprintf('\n initial kappa=...%+3.6f\n initial theta=...%+3.6f\n initial sigma = %+3.6f\n',kappa, theta, sigma);
if ~isfield(Model, 'MatlabDisp'), Model.MatlabDisp = 'off';
options = optimset('LargeScale', 'off', 'MaxIter', 300, 'MaxFunEvals',300, 'Display', Model.MatlabDisp, 'TolFun', 1e-4, 'TolX', 1e-4, 'TolCon', 1e-4); 
if ~isfield(Model, 'Method'), Model.Method = 'besseli'; 
if strcmp(Model.Method, 'ncx2pdf')
    [Params, Fval, Exitflag] =  fminsearch(@(Params) CIRobjective2(Params, Model),InitialParams, options);   
    [Params, Fval, Exitflag] =  fminsearch(@(Params) CIRobjective1(Params, Model),InitialParams, options);   
Results.Params = Params;
Results.Fval = -Fval/Nobs;
Results.Exitflag = Exitflag;

if strcmp(Model.Disp, 'y')
    fprintf('\n kappa = %+3.6f\n theta    = %+3.6f\n sigma = %+3.6f\n',...
        Params(1), Params(2), Params(3));
    fprintf(' log-likelihood = %+3.6f\n', -Fval/Nobs);



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.