these are the parameters that I use for the sort rateR0=0.05;sigma=0.1;alpha=0.6;bita=0.06;nSim=1000;N=360;%#of time stepsR=repmat(R0,nSim,1);%Simulation of the short ratesfor i=1:N R(:,i+1)=R(:,i)+alpha*(bita-R(:,i))*(1/N)+sigma*sqrt(R(:,i))*sqrt(1/N).*mvnrnd(0,1,nSim); endSo there are the short rates that I get, now I want to estimate the term structure using Rj=E((product(1+ri))^(1/j))y=cumprod(1+R');for i=1:N+1 y(i,=y(i,.^(1/i);endand then I use the closed form solution to the yield curver0=0.05;sigma=0.1;alpha=0.6;bita=0.06;P0=1;T=1;t=0;N=30;gamma=sqrt(alpha^2+2*sigma^2);for T=1:N B(T)=2*(exp(gamma*(T-t))-1)/((gamma+alpha)*(exp(gamma*(T-t))-1)+2*gamma); A(T)=(2*gamma*(exp((alpha+gamma)*(T-t)/2))/((gamma+alpha)*(exp(gamma*(T-t))-1)+2*gamma))^(2*alpha*bita/sigma^2); P(T)=A(T)*exp(-B(T)*r0); R(T)=-log(P(T))/(T-t); if T>1 f(T)=-log(P(T)/P(T-1)); else f(T)=-log(P(T)); end endbut the theoretical term structure is very different from the empirical...any sggestions????? Am I using the wrong formulas? is it a stupid Matlab bug? Please help

Somebody, please help.....

are the parameters sigma, alpha, bita correctly estimated from real rate?alpha seems too big...Bill

Hi,The parameters are just dummy parameters....

If I got your point, I don´t have a good news to you. CIR models are not able to match empirical Term Structure. To do that you need use HJM framework, try for instance, HO & Lee or Hull & WHite

In the CIR framework, there is a closed formula for the term structure.In order to "test my code" I am trying to estimate the "model" term structure using the short rates, generated by the CIR PDE, and curves are very different

As I told you, CIR Models doesnt match empirical Spot Rate Curve

Could you elaborate on that?Thank you so much!

Hi there, the code below seems to be working. There were some inconsistencies with your discounting intervals. Let me know if you need any clarifications.Kyriakosclear;clc;R0 = 0.05;sigma = 0.1;alpha = 0.6;bita = 0.06;nSim = 10000;N = 360;dt = 1/N;R = R0*ones(N+1,nSim);%Simulatedfor i=1:N W = randn(1,nSim/2); R(i+1,: ) = R(i,: ) + dt*alpha*(bita-R(i,: ))... +sigma*sqrt(dt*R(i,: )).*[W, -W];endy = 1./cumprod(1+dt*R);y0 = mean(y')';T0 = dt*(1:N+1)';P0 = -log(y0)./T0;plot(T0,P0); hold on;%Closed formT=1;N=30;gamma=sqrt(alpha^2+2*sigma^2);for T=1:N TT = T/N; B(T)=2*(exp(gamma*TT)-1)/((gamma+alpha)*(exp(gamma*TT)-1)+2*gamma); A(T)=(2*gamma*(exp((alpha+gamma)*TT/2))/((gamma+alpha)*(exp(gamma*TT)-1)+2*gamma))^(2*alpha*bita/sigma^2); P(T)=A(T)*exp(-B(T)*R0); Y(T)=-log(P(T))/TT; if T>1 f(T)=-log(P(T)/P(T-1)); else f(T)=-log(P(T)); endendplot((1:N)./N, Y,'r*');axis([0 1 -Inf Inf]);grid on;