Maximum Likelihood Estimation with Kalman filter using fminsearch

I have written a function to estimate a state-space model (Gaussian Affine Term Structure Model) via Maximum Likelihood Estimation with a Kalman filter. I am using fminsearch to optimize the log-likelihood, but for some unknown reasons it does not optimize anything and does not return any error. The returned vector of parameters Omega is my first guess...

function [Omega]=GATSM

      ns=3; %Number variable in the state vector
      maturities=[1,2,4,8,12,20,28,40,80,120]; %Maturities of the bond yields in quarters
      nyields=size(maturities,2);
      ytm=zeros(10,30);
      T = size(ytm,2);

      % Parameterization
      % Size of the matrices
      muP=NaN(ns,1);
      thetaP = NaN(ns,ns);
      muQ=NaN(ns,1);
      thetaQ = NaN(ns,ns);
      sigma_epsilon=NaN(ns,ns);
      delta0=NaN(1,1);
      delta1=NaN(ns,1);

      % Kalman Filter %
      % STEP 0: Setup
      Omega(1,:)=[0.5,0.25,0.36,0.42,0.5,0.2,0.3,0.5,2.5,3.1,0.6,0.2,0.4,0.88,0.69,0.41,0.2,0.12,0.10,0.5,0.6,0.3,0.2];

for i=1:10    
    % First guess and parameters restrictions
    sigma_epsilon=[Omega(i,17),0,0;Omega(i,18),Omega(i,19),0;Omega(i,20),Omega(i,21),Omega(i,22)];
    muP=[Omega(i,2);Omega(i,3);Omega(i,4)];
    thetaP = [Omega(i,8),Omega(i,9),Omega(i,10);Omega(i,11),Omega(i,12),Omega(i,13);Omega(i,14),Omega(i,15),Omega(i,16)];
    muQ=[Omega(i,1);0;0];
    thetaQ = diag([Omega(i,5) Omega(i,6) Omega(i,7)]);
    varsigma=Omega(i,23);
    delta0=0;
    delta1=[1;1;1];
    sigma_varsigma=varsigma*eye(nyields);

      % Compute A and B (Gaussian Affine Term Structure Model)
      A = NaN(nyields,1);
      B = NaN(ns,nyields);
      An=NaN(120,1);
      Bn=NaN(ns,120);
      An(1)=-delta0;
      Bn(:,1)=-delta1;
      for j = 2:1:120 % 120Q being the maturity of the longer bond
          Bn(:,j)= -(transpose(Bn(:,j-1))*thetaQ-transpose(delta1))/j;
          An(j) = -(An(j-1) + transpose(Bn(:,j-1))*muQ-delta0+0.5*transpose(Bn(:,j-1))*sigma_epsilon*transpose(sigma_epsilon)*Bn(:,j-1))/j;
      end
      A=[An(1);An(2);An(4);An(8);An(12);An(20);An(28);An(40);An(80);An(120)];
      B=[Bn(:,1),Bn(:,2),Bn(:,4),Bn(:,8),Bn(:,12),Bn(:,20),Bn(:,28),Bn(:,40),Bn(:,80),Bn(:,120)];

      % STEP 1: Initialization
      X1 = NaN(ns,T+1); % E_t-1[X_t]; t = 1,...,T
      X2 = NaN(ns,T+1); % E_t-1[X_t-1]; t = 1,...,T; only meanningful after t>=2
      P1 = NaN(ns,ns,T+1); % COV_t-1[X_t]; t = 1,...,T
      P2 = NaN(ns,ns,T+1); % COV_t-1[X_t-1]; t = 1,...,T; only meanningful after t>=2
      X2(:,1) = zeros(ns,1); % Come from Wu & Xia code
      P2 (:,:,1) = 100*eye(ns);
      loglikvec=NaN(T,1);

      % Filtering
      for t = 2:1:T+1
          % STEP 2: Predict
          X1(:,t) = muP + thetaP*X2(:,t-1);
          P1(:,:,t-1) = thetaP*P2(:,:,t-1)*thetaP' + sigma_epsilon*sigma_epsilon';
          ytm_hat = A + B'*X1(:,t);
          err = ytm(:,t-1) - ytm_hat;
          F = B'*P1(:,:,t-1)*B + sigma_varsigma*sigma_varsigma';

          % STEP 3: Update
          K = P1(:,:,t-1)*B*inv(F); % Kalman Gain
          X2(:,t) = X1(:,t) + K*err;
          P2(:,:,t) = P1(:,:,t-1)-K*B'*P1(:,:,t-1);

          % STEP 5: Calculate the likelihood function
          loglikvec(t) = -0.5*(nyields*log(2*pi) + log(det(F)) + err'*inv(F)*err);
      end
      llf = -sum(loglikvec);
      Omeo=Omega(i,:);
      fun=@(Omeo)llf;
      Omega(i+1,:) = fminsearch(fun,Omeo);
  end

      implied_X = X2(:,2:end);
      implied_Y = repmat(A,1,T) + B'*implied_X;

end

Thanks in advance for your precious help!

Answer 1

[Matt J]

Your function

 fun=@(Omeo)llf;

has no dependence on Omeo. It will return the constant vector llf all the time.

As a side note, fminsearch is wrong for this problem. You have way too many unknowns.

[Matt J]

Etienne Vaccaro-Grange commented:

Thanks for your prompt reply Matt.

My function actually depends on Omeo but indirectly (through other variables).

How can I specify my function knowing that it is calculated in a for loop with different steps?

[Matt J]

How can I specify my function knowing that it is calculated in a for loop with different steps?

You need to write a function file or a subfunction that executes those steps in its own workspace (not in the workspace of your main function).

I think it will also help you to study the following two examples and to understand why fun() always returns 4 in this first example,

>> omega=1; x=omega+3; fun=@(omega) x;

>> fun(1)

ans =

       4

>> fun(2)

ans =

       4

>> fun(3)

ans =

       4

but has different behavior in this 2nd example,

>> b=3; fun=@(omega) omega+b;

>> fun(1)

ans =

       4

>> fun(2)

ans =

       5

>> fun(3)

ans =

       6

[Matt J]

Etienne Vaccaro-Grange commented:

Thank you Matt, now it works.

As you advised, I created another function in another script calling the function to optimize.

I also changed fminsearch for fminunc