function [M,Ma,N,Na] = swapsystem(this,ialt,exi,endi,last)
% swapsystem  [Not a public function] Model solution matrices with some of the variables exogenised and shocks endogenised.
%
% Backend IRIS function.
% No help provided.

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

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

ny = sum(this.nametype == 1);
nx = size(this.solution{1},1);
nb = size(this.solution{1},2);
nf = nx - nb;
ne = sum(this.nametype == 3);

% Current-dated variables in the original state vector.
xcurri = imag(this.solutionid{2}) == 0;
nxcurr = sum(xcurri);
fcurri = xcurri(1:nf);
bcurri = xcurri(nf+1:end);

% Constant.
Mxc = zeros(nxcurr*last,1);
Myc = zeros(ny*last,1);
mac = zeros(nb,1);
% Multipliers on initial condition.
Mx0 = zeros(nxcurr*last,nb);
My0 = zeros(ny*last,nb);
ma0 = eye(nb);
% Multipliers on unexpected shocks.
Mxu = zeros(nxcurr*last,ne*last);
Myu = zeros(ny*last,ne*last);
mau = zeros(nb,ne*last);
% Multipliers on expected shocks.
Mxe = zeros(nxcurr*last,ne*last);
Mye = zeros(ny*last,ne*last);
mae = zeros(nb,ne*last);

% System matrices.
Tf = this.solution{1}(1:nf,:,ialt);
Ta = this.solution{1}(nf+1:end,:,ialt);
Ru = this.solution{2}(:,1:ne,ialt);
Re = this.solution{2}(:,1:ne*last,ialt);
Kf = this.solution{3}(1:nf,ialt);
Ka = this.solution{3}(nf+1:end,ialt);
Z = this.solution{4}(:,:,ialt);
H = this.solution{5}(:,:,ialt);
D = this.solution{6}(:,ialt);
U = this.solution{7}(:,:,ialt);
Ucurr = U(bcurri,:);

for t = 1 : last
    % Constant.
    mfc = Tf*mac + Kf;
    mac = Ta*mac + Ka;
    Mxc((t-1)*nxcurr+(1:nxcurr),1) = [mfc(fcurri,:);Ucurr*mac];
    Myc((t-1)*ny+(1:ny),1) = Z*mac + D;
    % Initial condition.
    mf0 = Tf*ma0;
    ma0 = Ta*ma0;
    Mx0((t-1)*nxcurr+(1:nxcurr),:) = [mf0(fcurri,:);Ucurr*ma0];
    My0((t-1)*ny+(1:ny),:) = Z*ma0;
    % Unexpected.
    mfu = Tf*mau;
    mau = Ta*mau;
    mfu(:,(t-1)*ne+(1:ne)) = mfu(:,(t-1)*ne+(1:ne)) + Ru(1:nf,:);
    mau(:,(t-1)*ne+(1:ne)) = mau(:,(t-1)*ne+(1:ne)) + Ru(nf+1:end,:);
    myu = Z*mau;
    myu(:,(t-1)*ne+(1:ne)) = myu(:,(t-1)*ne+(1:ne)) + H;
    Mxu((t-1)*nxcurr+(1:nxcurr),:) = [mfu(fcurri,:);Ucurr*mau];
    Myu((t-1)*ny+(1:ny),:) = myu;
    % Expected.
    mfe = Tf*mae;
    mae = Ta*mae;
    mfe(:,(t-1)*ne+1:end) = mfe(:,(t-1)*ne+1:end) + Re(1:nf,:);
    mae(:,(t-1)*ne+1:end) = mae(:,(t-1)*ne+1:end) + Re(nf+1:end,:);
    Re(:,end-ne+1:end) = [];
    mye = Z*mae;
    mye(:,(t-1)*ne+(1:ne)) = mye(:,(t-1)*ne+(1:ne)) + H;
    Mxe((t-1)*nxcurr+(1:nxcurr),:) = [mfe(fcurri,:);Ucurr*mae];
    Mye((t-1)*ny+(1:ny),:) = mye;
end

% Original system I*[Y;X] = M1*[U;E].
M = [Myc,My0,Myu,Mye;Mxc,Mx0,Mxu,Mxe];
Ma = [mac,ma0,mau,mae];

% When computing MSE matrices, we treat expected shocks as unexpected.
if nargout > 2
    N = [Myc,My0,Myu,Myu;Mxc,Mx0,Mxu,Mxu];
    Na = [mac,ma0,mau,mau];
end

if any(exi) || any(endi)
    % Add the alpha vector at t=last so that it is easy to retrieve the
    % initial condition for simulating the model after t=last.
    M = [Ma;M];
    exi = [false(1,nb),exi];
    % Swap the endogenised and exogenised columns in I and A matrices.
    I = eye(size(M,1));
    I1 = I(:,~exi);
    I2 = I(:,exi);
    M1 = M(:,~endi);
    M2 = M(:,endi);
    M = [I1,-M2]\[M1,-I2];
    Ma = M(1:nb,:);
    M = M(nb+1:end,:);
    if nargout > 2
        N = [Na;N];
        N1 = N(:,~endi);
        N2 = N(:,endi);
        N = [I1,-N2]\[N1,-I2];
        Na = N(1:nb,:);
        N = N(nb+1:end,:);
    end
end

end