wolffd@0: function [ypred, ll, mse] = eval_AR_perf(coef, C, y, model)
wolffd@0: % Evaluate the performance of an AR model.
wolffd@0: % 
wolffd@0: % Inputs
wolffd@0: % coef(:,:,k,m) - coef. matrix to use for k steps back, model m
wolffd@0: % C(:,:,m)      - cov. matrix for model m
wolffd@0: % y(:,t)        - observation at time t
wolffd@0: % model(t)      - which model to use at time t (defaults to 1 if not specified)
wolffd@0: %
wolffd@0: % Outputs
wolffd@0: % ypred(:,t)    - the predicted value of y at t based on the evidence thru t-1.
wolffd@0: % ll            - log likelihood
wolffd@0: % mse           - mean squared error = sum_t d_t . d_t, where d_t = pred(y_t) - y(t)
wolffd@0: 
wolffd@0: [s T] = size(y);
wolffd@0: k = size(coef, 3);
wolffd@0: M = size(coef, 4);
wolffd@0: 
wolffd@0: if nargin<4, model = ones(1, T); end
wolffd@0: 
wolffd@0: ypred = zeros(s, T);
wolffd@0: ypred(:, 1:k) = y(:, 1:k);
wolffd@0: mse = 0;
wolffd@0: ll = 0;
wolffd@0: for j=1:M
wolffd@0:   c(j) = log(normal_coef(C(:,:,j)));
wolffd@0:   invC(:,:,j) = inv(C(:,:,j));
wolffd@0: end
wolffd@0: coef = reshape(coef, [s s*k M]);
wolffd@0: 
wolffd@0: for t=k+1:T
wolffd@0:   m = model(t-k);
wolffd@0:   past = y(:,t-1:-1:t-k);
wolffd@0:   ypred(:,t) = coef(:, :, m) * past(:);
wolffd@0:   d = ypred(:,t) - y(:,t);
wolffd@0:   mse = mse + d' * d;
wolffd@0:   ll = ll + c(m) - 0.5*(d' * invC(:,:,m) * d);
wolffd@0: end
wolffd@0: mse = mse / (T-k+1);
wolffd@0: