function [l,omega,pstar,Aobs,exitflagFP,grate] = LogLikelihood(params,Pobs,model,interp,options,purge)
% LOGLIKELIHOOD Calculates the log-likelihood for given price observations and given parameters

%% Initialization
if nargin<6, purge = false; end
persistent cx x
if purge
  cx = [];
  x  = [];
end
if ~isempty(cx), interp.cx = cx; end
if ~isempty(x),  interp.x  = x;  end

model.params(1:4)                    = params(end-3:end);
par                                  = num2cell(model.params);
g                                    = params(1:end-4);
[a, b, delta]                        = par{1:3};
[interp,exitflagFP,pstar]            = SolveStorageDL(model,interp,options);
cx                                   = interp.cx;
A                                    = interp.s;
x                                    = interp.x;
T                                    = length(Pobs);
grate                                = model.basis*g;
Psto                                 = Pobs.*exp(-grate);
demand                               = @(p) (p-a)/b;

astar = demand(pstar);
if astar <= interp.s(1)
  warning('Availability cutoff below lower bound, astar= %0.3g%',astar)
end

%% Find availabilities corresponding to observed prices
invPriceFunction = interp1(x(:,2),A,options.InterpMethod,'pp');
Aobs             = demand(Psto);
infPstar         = Psto<pstar;
Aobs(infPstar)   = max(ppval(invPriceFunction,Psto(infPstar)),Aobs(infPstar));

%% Residuals
Sobs       = max(Aobs-demand(Psto),0);
omega      = NaN(T,1);
omega(2:T) = Aobs(2:T) - (1-delta)*Sobs(1:T-1);
OmegaExceedBounds = abs(omega)>5;

%% Jacobian
%{
There are three methods to calculate the Jacobian:
 1. Calculate the derivative of the inverse price function by finite
    differences.
 2. Use the fact that (f^{-1})'(f(x))=1/f'(x) and calculate the inverse of the
    derivative of the price function.
 3. Calculate the derivative of the inverse price function analytically.
In all cases, we can correct for the fact that above pstar the derivative of
the inverse price function is equal to 1/b and overall it never exceeds 1/b.

The methods 1 and 3 are equivalent, except that 3 is faster and more precise.
The methods 2 and 3 are equivalent for high levels of precision of the spline
approximation. We opt for method 3, which should be the fastest method, since
there is one less operation.
%}
% 1. Finite-difference of the inverse price function
% J              = min(diag(numjac(@(P) ppval(invPriceFunction,P),Pobs)),1/b);
% J(Pobs>=pstar) = 1/b;

% 2. Inverse of the derivative of the price function
% [breaks,coefs,l,order,d] = unmkpp(interp.cx{2});
% dPriceFunction           = mkpp(breaks,repmat(order-1:-1:1,d*l,1).*coefs(:,1:order-1),d);
% J                        = ones(size(Pobs))/b;
% J(Pobs<pstar)            = min(1./dPriceFunction(Aobs(Pobs<pstar)),1/b);

% 3. Differentiate the inverse price function
[breaks,coefs,l,order,d] = unmkpp(invPriceFunction);
dinvPriceFunction        = mkpp(breaks,repmat(order-1:-1:1,d*l,1).*coefs(:,1:order-1),d);
J                        = ones(size(Psto))/b;
J(infPstar)              = min(ppval(dinvPriceFunction,Psto(infPstar)),1/b);

J                        = exp(-grate).*J;


%% Log-likelihood
switch options.estimator
  case 'UML'
    %% Probability of the first observation
    Prob1                = UnconditionalProbability(Aobs(1),J(1),model,interp);

    %% Log-likelihood
    l = [log(Prob1)
         -log(model.TruncCorrection) - ...
         0.5*(log(2*pi) + ...
              (~OmegaExceedBounds(2:T)).*omega(2:T).^2+OmegaExceedBounds(2:T)*realmax) + ...
         log(abs(J(2:T)))];
  case 'CML'
    %% Log-likelihood
    l = -log(model.TruncCorrection) - ...
        0.5*(log(2*pi) + ...
             (~OmegaExceedBounds(2:T)).*omega(2:T).^2+OmegaExceedBounds(2:T)*realmax) + ...
        log(abs(J(2:T)));
end
