function llfn = llfn(bigthet);
% Uses the Kalman filter to evaluate the negative log likelihood
% function for the model with sticky prices. The parameters
% are transformed so that they satisfy theoretical restrictions.
%
% THIS PROGRAM WAS WRITTEN FOR MATLAB BY
%
%   PETER N. IRELAND
%   BOSTON COLLEGE
%   DEPARTMENT OF ECONOMICS
%   140 COMMONWEALTH AVENUE
%   CHESTNUT HILL, MA 02467
%   irelandp@bc.edu
%
%  FINANCIAL SUPPORT FROM THE NATIONAL SCIENCE FOUNDATION UNDER
%    GRANT NO. SES-0213461 IS GRATEFULLY ACKNOWLEDGED.
%
%  COPYRIGHT (c) 2002 BY PETER N. IRELAND.  REDISTRIBUTION IS
%    PERMITTED FOR EDUCATIONAL AND RESEARCH PURPOSES, SO LONG AS
%    NO CHANGES ARE MADE.  ALL COPIES MUST BE PROVIDED FREE OF
%    CHARGE AND MUST INCLUDE THIS COPYRIGHT NOTICE.
%
%  Modified by Oke Roehe (2011)

% define variables and parameters
load Slovenia_centered2;

  global ct it mt pt rt bigt

  bigthet = real(bigthet);

  betatr = bigthet(1);
  gammatr = bigthet(2);
  etatr = 1.5;
  thetatr = 5;
  alphatr = 0.29;
  gtr = 0;
  deltatr = sqrt(0.025/0.975);

  phiptr = bigthet(3);
  phiktr = bigthet(4);

  mutr = bigthet(5);
  omegartr =1;
  omegamtr = bigthet(6);
  omegaptr = bigthet(7);
  omegaytr = bigthet(8);
    
  etr = bigthet(9);
  ztr = bigthet(10);

  rhoatr = bigthet(11);
  rhoetr = bigthet(12);
  rhoxtr = bigthet(13);
  rhoztr = bigthet(14);
  rhovtr = bigthet(15);

  sigatr = bigthet(16);
  sigetr = bigthet(17);
  sigxtr = bigthet(18);
  sigztr = bigthet(19);
  sigvtr = bigthet(20);
  
  %untransform parameters

  beta = (100*betatr)^2/(1+(100*betatr)^2);
  gamma = gammatr^2/(1+gammatr^2);% du to random process
  eta = abs(etatr);
  theta = 1 + abs(thetatr);
  alpha = abs(alphatr);
  g = 1 + abs(gtr);
  delta = deltatr^2/(1+deltatr^2);
  
  phip = abs(phiptr)*100;
  phik = abs(phiktr)*100;

  mu = 1 + abs(mutr);
  omegar = omegartr;
  omegam = omegamtr;
  omegap = omegaptr;
  omegay = omegaytr;

  e = abs(etr);
  z = 10000*abs(ztr);

  rhoa = (10*rhoatr)^2/(1+(10*rhoatr)^2);
  rhoe = (10*rhoetr)^2/(1+(10*rhoetr)^2);
  rhox = (10*rhoxtr)^2/(1+(10*rhoxtr)^2);
  rhoz = (10*rhoztr)^2/(1+(10*rhoztr)^2);
  rhov = (10*rhovtr)^2/(1+(10*rhovtr)^2);

  siga = abs(sigatr);
  sige = abs(sigetr);
  sigx = abs(sigxtr);
  sigz = abs(sigztr);
  sigv = abs(sigvtr);

% compute steady-state values

  piss = mu/g;

  rss = mu/beta;

  qss = g/beta - 1 + delta;

  l1ss = 1/(theta/(theta-1)-(g-1+delta)*(alpha/qss));

  l2ss = 1/(1+e*(rss/(rss-1))^(gamma-1));

  l3ss = (1-alpha)*z*((theta-1)/theta)^(1/(1-alpha)) ...
           *(alpha/qss)^(alpha/(1-alpha));

  lambdass = (eta+(1-alpha)*l1ss*l2ss)/l3ss;

  xiss = ((theta-1)/theta)*lambdass;

  css = (1/lambdass)/(1+e*(rss/(rss-1))^(gamma-1));

  mss = e*((rss/(rss-1))^gamma)*css;

  yss = css/(1-(g-1+delta)*((theta-1)/theta)*(alpha/qss));

  kss = ((theta-1)/theta)*alpha*yss/qss;

  iss = (g-1+delta)*kss;

  hss = (1/z)*((yss/(kss^alpha))^(1/(1-alpha)));

  wss = (1-alpha)*((theta-1)/theta)*yss/hss;

  dss = yss - wss*hss - qss*kss;

  nss = yss/hss;

% form matrices A, B, and C

  biga = zeros(10,10);
  bigb = zeros(10,5);
  bigc = zeros(10,5);

  biga(1,1) = yss;
  biga(1,2) = -css;
  biga(1,3) = -iss;

  biga(2,2) = rss*(1+(gamma-1)*lambdass*css);
  biga(2,6) = (gamma-1)*(rss-1)*lambdass*mss;
  bigb(2,2) = -(gamma-1)*(rss-1)*lambdass*mss;
  bigb(2,3) = (gamma-1)*(rss-1)*lambdass*mss;
  bigb(2,4) = -gamma*rss;
  bigc(2,1) = gamma*rss;
  bigc(2,2) = -(rss-1)*lambdass*mss;

  biga(3,4) = lambdass*wss*hss;
  biga(3,7) = -eta;
  bigb(3,4) = eta;

  biga(4,2) = rss - 1;
  biga(4,6) = -(rss-1);
  biga(4,10) = -gamma;
  bigb(4,2) = rss - 1;
  bigb(4,3) = -(rss-1);
  bigc(4,2) = -(rss-1);

  biga(5,1) = yss;
  biga(5,4) = -wss*hss;
  biga(5,7) = -wss*hss;
  biga(5,8) = -qss*kss;
  biga(5,9) = -dss;
  bigb(5,1) = qss*kss;

  biga(6,1) = 1;
  biga(6,4) = -(1-alpha);
  bigb(6,1) = alpha;
  bigc(6,4) = 1 - alpha;

  biga(7,1) = 1;
  biga(7,4) = -1;
  biga(7,7) = -1;
  bigb(7,4) = 1;
  bigb(7,5) = -1;

  biga(8,1) = 1;
  biga(8,8) = -1;
  bigb(8,1) = 1;
  bigb(8,4) = 1;
  bigb(8,5) = -1;

  biga(9,1) = omegay;
  biga(9,6) = omegam;
  biga(9,10) = -omegar;
  bigb(9,3) = -omegap;
  bigc(9,5) = -1;

  biga(10,1) = 1;
  biga(10,4) = -1;
  biga(10,5) = -1;

% form matrices D, F, G, H, and J

  bigd = zeros(5,5);
  bigf = zeros(5,10);
  bigg = zeros(5,5);
  bigh = zeros(5,10);
  bigj = zeros(5,5);

  bigd(1,1) = g*kss;
  bigg(1,1) = (1-delta)*kss;
  bigh(1,3) = iss;
  bigj(1,3) = iss;

  bigd(2,3) = 1;
  bigd(2,4) = -1;
  bigg(2,4) = -1;
  bigh(2,10) = 1;

  bigd(3,1) = phik*(beta*(1-delta)/g-(1+beta));
  bigd(3,4) = g;
  bigf(3,3) = beta*phik*(1-(1-delta)/g);
  bigf(3,8) = beta*qss;
  bigg(3,1) = -phik;
  bigg(3,4) = g;
  bigj(3,3) = -g-beta*(phik*(1-(1-delta)/g)-(1-delta))*rhox;

  bigd(4,3) = beta*phip;
  bigg(4,3) = phip;
  bigg(4,4) = theta - 1;
  bigg(4,5) = -(theta-1);

  bigd(5,2) = 1;
  bigg(5,2) = 1;
  bigg(5,3) = -1;
  bigh(5,6) = 1;

% form matrix P

  bigp = zeros(5,5);

  bigp(1,1) = rhoa;
  bigp(2,2) = rhoe;
  bigp(3,3) = rhox;
  bigp(4,4) = rhoz;
  bigp(5,5) = rhov;

% form matrices K and L

  bigkl1 = inv(bigd+bigf*inv(biga)*bigb);
  bigkl2 = bigg + bigh*inv(biga)*bigb;
  bigkl3 = bigj + bigh*inv(biga)*bigc - bigf*inv(biga)*bigc*bigp;

  bigk = bigkl1*bigkl2;

  bigl = bigkl1*bigkl3;

  bigl1 = bigl(1:2,:);
  bigl2 = bigl(3:5,:);

% form matrices M and N

  [kvec,keig] = eig(bigk);

  keig2 = diag(keig);

  [keig3,kord] = sort(abs(keig2));

  kviol = 0;

  if keig3(2) > 1
    kviol = 1;
  end

  if keig3(3) < 1
    kviol = 1;
  end

  keig4 = keig2(kord);

  bign = diag(keig4);

  bign1 = diag(keig4(1:2));
  bign2 = diag(keig4(3:5));

  bigminv = kvec(:,kord);

  bigm = inv(bigminv);

  bigm11 = bigm(1:2,1:2);
  bigm12 = bigm(1:2,3:5);
  bigm21 = bigm(3:5,1:2);
  bigm22 = bigm(3:5,3:5);

% form matrices Q and R

  bigq1 = bigm11*bigl1 + bigm12*bigl2;
  bigq2 = bigm21*bigl1 + bigm22*bigl2;

  vecr = inv(eye(15)-kron(bigp,inv(bign2)))*bigq2(:);

  bigr = reshape(vecr,3,5);

% form matrices S

  bigs1 = -inv(bigm22)*bigm21;

  bigs2 = -inv(bigm22)*inv(bign2)*bigr;

  bigs34 = bigm11 + bigm12*bigs1;

  bigs3 = inv(bigs34)*bign1*bigs34;

  bigs4 = inv(bigs34)*(bigq1+bign1*bigm12*bigs2-bigm12*bigs2*bigp);

  bigs5a = [ eye(2) ; bigs1 ]; 

  bigs5 = inv(biga)*bigb*bigs5a;

  bigs6a = [ zeros(2,5) ; bigs2 ];

  bigs6 = inv(biga)*bigc + inv(biga)*bigb*bigs6a;

% form matrices PI, W, and U

  bigpi = [ bigs3 bigs4 ; zeros(5,2) bigp ];

  bigw = [ zeros(2,5) ; eye(5) ];

  bigu = [ bigs5 bigs6 ; bigs1 bigs2 ];

  bigpi = real(bigpi);

  bigu = real(bigu);

% form matrices AX, BX, CX, VX, and BVBX

  bigax = bigpi;

  bigbx = bigw;

  bigcx = [ bigu(2,:) ; ...
            bigu(3,:) ; ...
            bigpi(2,:) ; ...
            bigu(11,:) ; ...
            bigu(10,:) ];

  bigvx = diag([siga^2 sige^2 sigx^2 sigz^2 sigv^2]);

  bigbvbx = bigbx*bigvx*bigbx';

% put data in deviation form

  trend = 1:bigt;

%  trend = trend + 82;

  cthat = ct - log(g)*trend' - log(css);  % Keep this structure to allow for linear detrended data (just in case) -> here: log(g)=0

  ithat = it - log(g)*trend' - log(iss);

  mthat = mt - log(g)*trend' - log(mss);

  pthat = pt - log(piss);

  rthat = rt - log(rss);

  dthat = [ cthat ithat mthat pthat rthat ];

% evaluate negative log likelihood

  st = zeros(7,1);

  bigsig1 = inv(eye(49)-kron(bigax,bigax))*bigbvbx(:);

  bigsigt = reshape(bigsig1,7,7);

  llfn = (5*bigt/2)*log(2*pi);

  for t = 1:bigt

    ut = dthat(t,:)' - bigcx*st;

    omegt = bigcx*bigsigt*bigcx';

    omeginvt = inv(omegt);

    llfn = llfn + (1/2)*(log(det(omegt))+ut'*omeginvt*ut);

    bigkt = bigax*bigsigt*bigcx'*omeginvt;

    st = bigax*st + bigkt*ut;

    bigsigt = bigbvbx + bigax*bigsigt*bigax' ...
                - bigax*bigsigt*bigcx'*omeginvt*bigcx*bigsigt*bigax';

  end

% penalize eigenvalue violations

  if kviol

    llfn = llfn + 1e12;

  end

  if abs(imag(llfn)) > 0

    llfn = real(llfn) + 1e12;

  end