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 '14 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$ – SchnitzelRaver Dec 23 '14 at 11:16
  • $\begingroup$ Here it is! Heston calibration in Matlab $\endgroup$ – Gabriele Pompa Dec 23 '14 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$ – Joshua Ulrich Dec 23 '14 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$ – SchnitzelRaver Dec 23 '14 at 16:37
  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, privacy policy and cookie policy

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