www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xregtwostage/evalpev.m

    function  PEVout = evalpev(X,TS,translocal)
%EVALPEV Evaluate prediction error variance for a twostage model
%
%  PEVout=EVALPEV(X,m);
%
%  This function is normally called via model/PEV
%
%  A function to compute the Prediction error Variance of a multistage
%  model. xg represents the values of the local and then global variables,
%  m is a Two Stage model object.

%  Copyright 2000-2007 The MathWorks, Inc. and Ford Global Technologies, Inc.


if nargin<3
    translocal= true;
end
L= TS.Local;
Nf= size(L,1);

R= var(TS);
if isempty(R)
    error(message('mbc:xregtwostage:InvalidState'));
end

% number of local variables
nl = nfactors(L);

if iscell(X)
    GridEval= true;
    % cell array used for pevgrid
    Xl= X{1};
    Xg= X{2};

    Ng= size(Xg,1);
    PEVout= cell(Ng,1);
else
    GridEval= false;
    Xl= X(:,1:nl) ;					% Get local variable values.
    Xg= X(:,nl+1:end);			    % Get global variable values.
    Ng= size(Xg,1);
    PEVout= zeros(Ng,1);
end

% global part of the jacobian
D= JacobGlobalVar(TS,Xg)*R;

% evaluate and get reconstruction parameters
[Y,Yg,Datum,Lparams] = eval(TS,X);

reEvalDG= ~allLinearRF(L);
if ~reEvalDG
    % inv(delg/delp) should be reevaluated here
    L= EvalDelG(L);
    dh= inv(delG(L));
end

if ~GridEval
    Lstore= cell(1,Ng);
end

for i= 1:Ng

    % Note that the last input here treats the local variable
    %  as relative to the datum
    % this only works for a single output
    L= datum(L,Datum(i,:));
    L= update(L,Lparams(i,:)',[]);

    if reEvalDG
        % inv(delg/delp) should be reevaluated here
        L= EvalDelG(L);
        dh= inv(delG(L));
    end

    % local jacobian
    if GridEval
        J = jacobian(L,Xl,0);
        
        % local sweep jacobian* inv(delg/delp)
        M= J*dh;

        % current sweep
        M= M*D((i-1)*Nf+1:Nf*i,:);
        
        PEVout{i}= sum(M.*M,2);
    else
        J = jacobian(L,Xl(i,:),0);
        Lstore{i}= J*dh;
    end
end

if GridEval
    PEVout= cat(1,PEVout{:});
else
    M= spblkdiag(Lstore{:})*D;
    PEVout= sum(M.*M,2);
end


if translocal && HasTransform(L)
    dy= yinvdiff(L,Y);
    PEVout= dy.^2.*PEVout;
end