Improving code execution: vectorization

I have these functions that calculate bond yields. In the code below, the most simple one is displayed, although there are similar function with more complicated code inside the loop. Execution time is of the essence as the functions go into larger optimization problems. The performance is highly dependent on the maturity input, which translates into the # of loops. Now, my question is whether it is somehow possible to rewrite this function - such as vectorization instead of looping - that will result in improvements of performance?

function [Yield] = Yield_shadowrate_FO(...
                                 alpha, beta, mu, phi, sigma, maturity, X)
  N = size(X, 1);
  s_mean = zeros(maturity-1, 1);
  s_var = zeros(maturity-1, 1);
  a = zeros(maturity-1, 1);
  FO = zeros(maturity-1, 1); 
  for k = 1:maturity
    if k == 1  
      short_rate = max(alpha+beta' * X,0);
    else
      i = (k-1);     
      X_mean = mu + (eye(N)-phi)^i*X - (phi\(eye(N)-phi)^i)*phi*mu;
      S = (eye(N*N)-kron(eye(N)-phi,eye(N)-phi)) \ ...
          (eye(N*N)-kron(eye(N)-phi,eye(N)-phi)^i)* ...
          reshape(sigma*sigma', N*N, 1);
      X_var = reshape(S, N, N);
      s_mean(i) = alpha + beta' * X_mean;
      s_var(i) = beta' * X_var * beta;
      a(i) = s_mean(i) / s_var(i)^0.5;
      FO(i)=s_mean(i)*0.5*erfc(-a(i)/sqrt(2))+s_var(i)^0.5*normpdf(a(i)); 
    end
  end
  Yield = 1/maturity*(short_rate+sum(FO,1)); 
end

[Geoff Hayes]

Anders - please comment on the performance of the above code and give us some indication of the size/dimension of each of the inputs alpha, beta, mu, phi, sigma, maturity, X. How often do you call the above function and how long does it take to execute it? Have you tried to profile the code?