function [this,data,B,COUNT] = myidentify(this,data,opt)
% myidentify  [Not a public function] Convert reduced-form VAR to structural VAR.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

%**************************************************************************

[ny,p,nalt] = sizeof(this); %#ok<ASGLU>
A = poly.var2poly(this.A);
Omega = this.Omega;

% Handle residuals.
isdata = nargin > 1 && nargout > 1 && ~isempty(data);

if isdata
    % Get data.
    [outputformat,range,y,e] = VAR.datarequest(this,data,Inf,opt);
    if size(e,3) == 1 && nalt > 1
        y = y(:,:,ones(1,nalt));
        e = e(:,:,ones(1,nalt));
    end
else
    y = [];
    e = [];
end

% Std dev of structural residuals requested by the user.
this.std(1,:) = opt.std(1,ones(1,nalt));

B = zeros(ny,ny,nalt);
COUNT = 1;
switch lower(opt.method)
    case 'chol'
        doreorder();
        for ialt = 1 : nalt
            B(:,:,ialt) = chol(Omega(:,:,ialt)).';
        end
        if opt.std ~= 1
            B = B / opt.std;
        end
        dobackorder();
        doconvertresid();
    case 'qr'
        doreorder();
        C = sum(A,3);
        for ialt = 1 : nalt
            B0 = transpose(chol(Omega(:,:,ialt)));
            if rank(C(:,:,1,ialt)) == ny
                Q = qr(transpose(C(:,:,1,ialt)\B0));
            else
                Q = qr(transpose(pinv(C(:,:,1,ialt))*B0));
            end
            B(:,:,ialt) = B0*Q;
        end
        if opt.std ~= 1
            B = B / opt.std;
        end
        dobackorder();
        doconvertresid();
    case 'svd'
        q = opt.rank;
        [B,e] = covfun.orthonorm(Omega,q,opt.std,e);
        % Recompute covariance matrix of reduced-form residuals if it is
        % reduced rank.
        if q < ny
            var = opt.std .^ 2;
            for ialt = 1 : nalt
                this.Omega(:,:,ialt) = ...
                    B(:,1:q,ialt)*B(:,1:q,ialt)'*var;
            end
            % Cannot produce structural residuals for reduced-rank cov matrix.
            data = [];
            isdata = false;
        else
            doconvertresid();
        end
    case 'householder'
        % Use Householder transformations to draw random SVARs. Test each SVAR
        % using `'test='` option to decide whether to keep it or discard.
        if nalt > 1
            utils.error('VAR', ...
                ['Cannot run SVAR() with ''method=search'' on ', ...
                'a multiple-param VAR object.']);
        end
        if isempty(opt.test)
            utils.error('VAR', ...
                ['Cannot run SVAR() with ''method=search'' and ', ...
                'empty ''test''.']);
        end
        if any(opt.ndraw <= 0)
            utils.warning('VAR', ...
                ['One of the ''ndraw'' options is zero, ', ...
                'returning an empty SVAR object.']);
        end
        [B,COUNT] = xxdraw(this,opt);
        nalt = size(B,3);
        if nalt > 0
            if opt.std ~= 1
                B = B / opt.std;
            end
            this = alter(this,nalt);
            if isdata
                % Expand the number of data sets to match the number of
                % newly created structural parameterisations.
                y = y(:,:,ones(1,nalt));
                e = e(:,:,ones(1,nalt));
            end
            doconvertresid();
        else
            isdata = false;
            this = alter(this,0);
            data = [];
        end
end

this.B(:,:,:) = B;
if isdata
    % Output data.
    names = get(this,'names');
    data = VAR.outputdata(this,outputformat,range,[y;e],[],names);
end

% Nested functions.

%**************************************************************************
    function doreorder()
        if ~isempty(opt.ordering)
            opt.ordering = opt.ordering(:)';
            if length(opt.ordering) ~= ny ...
                    || length(intersect(1:ny,opt.ordering)) ~= ny
                utils.error('VAR', ...
                    'Invalid ordering vector.');
            end
            A = A(opt.ordering,opt.ordering,:,:);
            Omega = Omega(opt.ordering,opt.ordering,:);
        end
    end
% doreorder().

%**************************************************************************
    function dobackorder()
        % Put variables (and residuals, if requested) back in order.
        if ~isempty(opt.ordering)
            [ans,backorder] = sort(opt.ordering); %#ok<NOANS,ASGLU>
            if opt.reorderresiduals
                B = B(backorder,backorder,:);
            else
                B = B(backorder,:,:);
            end
        end
    end
% dobackorder().

%**************************************************************************
    function doconvertresid()
        if isdata
            for iloop = 1 : nalt
                e(:,:,iloop) = B(:,:,iloop) \ e(:,:,iloop);
            end
        end
    end
% doconvertresid().

end

% Subfunctions.

%**************************************************************************
function [BB,COUNT] = xxdraw(this,opt)
%
% * Rubio-Ramirez J.F., D.Waggoner, T.Zha (2005) Markov-Switching Structural
% Vector Autoregressions: Theory and Application. FRB Atlanta 2005-27.
%
% * Berg T.O. (2010) Exploring the international transmission of U.S. stock
% price movements. Unpublished manuscript. Munich Personal RePEc Archive
% 23977, http://mpra.ub.uni-muenchen.de/23977.

test = opt.test;
A = poly.var2poly(this.A);
C = sum(A,3);
Ci = inv(C);
ny = size(A,1);

[h,isy] = mywoonvav(this,test);
P = covfun.orthonorm(this.Omega);
COUNT = 0;
maxfound = opt.ndraw;
maxiter = opt.maxiter;
BB = nan(ny,ny,0);
SS = nan(ny,ny,h,0);
YY = nan(ny,ny,0);

% Create command-window progress bar.
if opt.progress
    pbar = progressbar('IRIS VAR.SVAR progress');
end

nb = 0;
while COUNT < maxiter && nb < maxfound
    COUNT = COUNT + 1;
    % Candidate rotation. Note that we need to call `qr` with two
    % output arguments to get the unitary matrix `Q`.
    [Q,ans] = qr(randn(ny)); %#ok<NASGU,NOANS>
    B = P*Q;
    % Compute impulse responses T = 1 .. h.
    if h > 0
        S = timedom.var2vma(this.A,B,h);
    else
        S = zeros(ny,ny,0);
    end
    % Compute asymptotic cum responses.
    if isy
        Y = Ci*B; %#ok<MINV>
    end
    % Test impulse response.
    flag = false;
    dotest();
    nb = size(BB,3);
    if opt.progress
        update(pbar,max(COUNT/maxiter,nb/maxfound));
        % disp(max(count/maxcount,nb/maxfound));
    end
end

% Nested functions.

    function dotest()
        try
            flag = isequal(eval(test),true);
            if flag
                BB(:,:,end+1) = B;
                SS(:,:,:,end+1) = S; %#ok<SETNU>
                if isy
                    YY(:,:,end+1) = Y; %#ok<SETNU>
                end
            else
                % Test minus the structure.
                B = -B;
                S = -S;
                if isy
                    Y = -Y;
                end
                flag = isequal(eval(test),true);
                if flag
                    BB(:,:,end+1) = B;
                    SS(:,:,:,end+1) = S;
                    if isy
                        YY(:,:,end+1) = Y;
                    end
                end
            end
        catch err
            utils.error('VAR', ...
                ['Error evaluating the test string ''%s''.', ...
                '\nMatlab says ''%s'''], ...
                test,err.message);
        end
    end

end
% xxdraw().
