3
$\begingroup$

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.

Thanks!

Improve this question
$\endgroup$
7
  • $\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$
    – SchnitzelRaver
    Dec 23, 2014 at 11:16
  • $\begingroup$ Here it is! Heston calibration in Matlab $\endgroup$
    – Gabriele Pompa
    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$
    – Joshua Ulrich
    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$
    – SchnitzelRaver
    Dec 23, 2014 at 16:37

1 Answer 1

Reset to default
1
$\begingroup$
  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

Improve this answer
$\endgroup$

Your Answer

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

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