# Stochastic Volatility CIR estimation

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?

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

I just need to cross-check my results.

Thanks!

• 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. Commented Dec 23, 2014 at 10:52
• @Arrigo: Thanks, there's been a typo in the price process specification. Unlike for Heston, I need to estimate without using the option prices. Commented Dec 23, 2014 at 11:16
• Here it is! Heston calibration in Matlab Commented Dec 23, 2014 at 12:52
• Answers should contain enough information to actually answer the question. A simple link to another site should just be a comment. Commented Dec 23, 2014 at 14:08
• @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. Commented Dec 23, 2014 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';
end;
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);
end
if ~isfield(Model, 'MatlabDisp'), Model.MatlabDisp = 'off';
end;
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';
end;
if strcmp(Model.Method, 'ncx2pdf')
[Params, Fval, Exitflag] =  fminsearch(@(Params) CIRobjective2(Params, Model),InitialParams, options);
else
[Params, Fval, Exitflag] =  fminsearch(@(Params) CIRobjective1(Params, Model),InitialParams, options);
end
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);
end


end