![[FRB logo]](../../bitmaps/sm_head2.gif)
|
The Anderson-Moore Algorithm AIM Home Page |
|
Gary Anderson
[], []] describe a powerful method for solving linear saddle point models. The algorithm has proved useful in a wide array of applications including analyzing linear perfect foresight models, providing initial solutions and asymptotic constraints for nonlinear models. The algorithm solves linear problems with dozens of lags and leads and hundreds of equations in seconds. The technique works well for both symbolic algebra and numerical computation.
Although widely used at the Federal Reserve, few outside the central bank know about or have used the algorithm. This paper attempts to present the current algorithm in a more accessible format in the hope that economists outside the Federal Reserve may also find it useful. In addition, over the years there have been undocumented changes in approach that have improved the efficiency and reliability of algorithm. This paper describes the present state of development of this set of tools.
Anderson and Moore [[]] outlines a procedure that computes solutions for structural models of the form
The algorithm determines whether the model 1 has a unique solution, an infinity of solutions or no solutions at all.
The specification 1 is not restrictive.
One can handle inhomogeneous version of equation 1
by recasting the problem in terms of deviations from a steady state value or
by adding a new variable for each non-zero right hand side with an equation
guaranteeing the value always equals
the inhomogeneous value( 
and 
).
Saddle point problems combine initial conditions and asymptotic convergence to identify their solutions. The uniqueness of solutions to system 1 requires that the transition matrix characterizing the linear system have an appropriate number of explosive and stable eigenvalues[[]], and that the asymptotic linear constraints are linearly independent of explicit and implicit initial conditions[[]].
The solution methodology entails
Figure 1 presents a flow chart summarizing the algorithm. For a description of a parallel implementation see [[]] For a description of a continuous application see [[]].
The algorthim in this section performs elementary row
operations on the original model equations to produce
a normal form that facilitates construction of a state space
transition matrix.
If the leading square block of the linear system( 
) were non-singular,
one could compute the transition matrix from
![[AIM logo]](../../bitmaps/aimLogo1.0.gif)
"matlab/numericTransitionMatrix.m" 1 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function transMat=numericTransitionMatrix(transF)
[nr,nbc]=size(transF);
nc=nbc-nr;
%transMat = -transF(:,nc+1:nbc)\transF(:,1:nc);
transMat = -inv(transF(:,nc+1:nbc))*transF(:,1:nc);
if (nc > nr)
transMat=[zeros(nc-nr,nr),speye(nc-nr);transMat];
end
<>
blue
AN EXAMPLE:
Since, for most economic models,
The remaining equations, the auxiliary initial conditions,
provide important information for determining a unique trajectory for the
model.
Section 2.3 shows how
the auxiliary initial conditions provide a mechanism for reducing the
size of the invariant space calculation.
anewcolor
full(numericTransitionMatrix(sparse([1,2,3,4])))
full(numericTransitionMatrix(sparse([1,2,3,4])))
ans =
0 1.0000 0
0 0 1.0000
-0.2500 -0.5000 -0.7500
>>
>>
is typically singular, algorithm 1 identifies linear
combinations of the first 
equations that have leading
block non-singular. Inverting this non-singular right hand block
and multiplying the lagged matrices in the corresponding row
provides the autoregressive representation.
This routine checks the(dim) right most elements of a vector numbers(x)
to see if they are small. The routine
compares 
to 
and returns False if it
is larger and True otherwise.
"matlab/numericRightMostAllZeroQ.m" 2 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function allz = numericRightMostAllZeroQ(dim,x)
zeroTol=10^(-10);
allz=(norm(x(length(x)-dim+1:length(x))) < zeroTol);
<>
blue
AN EXAMPLE:
numericRightMostAllZeroQ(1,[1,2,3])
ans =
0
numericRightMostAllZeroQ(1,[1,2,0])
ans =
1
This routine checks the(dim) Left most elements of a vector numbers(x)
to see if they are small. The routine
compares 
to 
and returns False if it
is larger and True otherwise.
"matlab/numericLeftMostAllZeroQ.m" 3 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function allz = numericLeftMostAllZeroQ(dim,x)
zeroTol=10^(-10);
allz=(norm(x(1:dim)) < zeroTol);
<>
blue
AN EXAMPLE:
numericLeftMostAllZeroQ(1,[1,2,3])
ans =
0
numericLeftMostAllZeroQ(1,[0,1,2])
ans =
1
This routine implements lines 9-10 of the algorithm. It takes a set of auxiliary initial conditions and an input matrix(H) and
The routine should throw an exception if it identifies a matrix row that is all zeroes. Instead, the routine calle matlab's error function.
"matlab/numericShiftRightAndRecord.m" 4 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [ultimateAuxConds,shiftedHmat] = numericShiftRightAndRecord(auxConds,unshiftedHmat)
[hrows,hcols]=size(unshiftedHmat);
shiftedHmat=zeros(hrows,hcols);
rowsNow=0;
ultimateAuxConds=zeros(hcols-hrows,hcols-hrows);
zeroTol=10^(-10)*hrows;
for i = 1:hrows
if(norm(unshiftedHmat(i,:))<zeroTol)
error('zero row for hmat')
end;%if
while(numericRightMostAllZeroQ(hrows,unshiftedHmat(i,:)))
rowsNow=rowsNow+1;
ultimateAuxConds(rowsNow,:)=unshiftedHmat(i,1:hcols-hrows);
unshiftedHmat(i,:)=[zeros(1,hrows),unshiftedHmat(i,1:hcols-hrows)];
end; %end while
shiftedHmat(i,:)=unshiftedHmat(i,:);
end;%for,
ultimateAuxConds=[auxConds;ultimateAuxConds(1:rowsNow,:)];
<>
blue
AN EXAMPLE:
[ax,hf]=numericShiftRightAndRecord(sparse([1,2,7,9]),sparse([4,5,0,0,0,0;2,2,3,3,0,0]))
ax =
(1,1) 1
(2,1) 4
(4,1) 2
(1,2) 2
(2,2) 5
(4,2) 2
(1,3) 7
(3,3) 4
(4,3) 3
(1,4) 9
(3,4) 5
(4,4) 3
hf =
0 0 0 0 4 5
0 0 2 2 3 3
>>
This routine implements lines 
-10
of the algorithm.
It takes a set of auxiliary initial conditions and an input
matrix(H) and
The routine should throw an exception if it identifies a matrix row that is all zeroes. Instead, the routine calle matlab's error function.
"matlab/numericShiftLeftAndRecord.m" 5 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [ultimateAuxConds,shiftedHmat] = numericShiftLeftAndRecord(auxConds,unshiftedHmat)
[hrows,hcols]=size(unshiftedHmat);
shiftedHmat=zeros(hrows,hcols);
rowsNow=0;
ultimateAuxConds=zeros(hcols-hrows,hcols-hrows);
zeroTol=10^(-10)*hrows;
for i = 1:hrows
if(norm(unshiftedHmat(i,:))<zeroTol)
error('zero row for hmat')
end;%if
while(numericLeftMostAllZeroQ(hrows,unshiftedHmat(i,:)))
rowsNow=rowsNow+1;
ultimateAuxConds(rowsNow,:)=unshiftedHmat(i,hrows+1:hcols);
unshiftedHmat(i,:)=[unshiftedHmat(i,hrows+1:hcols),zeros(1,hrows)];
end; %end while
shiftedHmat(i,:)=unshiftedHmat(i,:);
end;%for,
ultimateAuxConds=[auxConds;ultimateAuxConds(1:rowsNow,:)];
<>
blue
AN EXAMPLE:
[ax,hf]=numericShiftLeftAndRecord(sparse([1,2,7,9]),sparse([0,0,4,5,0,0;0,0,0,0,3,3]))
ax =
(1,1) 1
(2,1) 4
(4,1) 3
(1,2) 2
(2,2) 5
(4,2) 3
(1,3) 7
(3,3) 3
(1,4) 9
(3,4) 3
hf =
4 5 0 0 0 0
3 3 0 0 0 0
>>
>>
This routine uses QR Decomposition with column pivoting to construct a
matrix that zeroes the rows of the input matrix.
The routine extends the matrix to a basis for the n dimensional
space by computing the orthogonal complement of 
.
"matlab/numericComputeAnnihilator.m" 6 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function annmat=numericComputeAnnihilator(amat,hmat)
zeroTol=10^(-10);
[arows,acols]=size(amat);
%p=arows:-1:1;
if(norm(reshape(amat,1,arows*acols))<zeroTol)
annmat=hmat;
else
if(issparse(amat))
[annmat,r]=qr(amat,hmat);
else
[annmat,r,e]=qr(amat);
annmat=(annmat')*hmat;
end%end if
%annmat=annmat(:,p);%want to place zero rows at top of matrix
%annmat=annmat';
end;%end if
<>
blue
AN EXAMPLE:
numericComputeAnnihilator(zeros(3),zeros(3))
ans =
0 0 0
0 0 0
0 0 0
>>
>>
>>
numericComputeAnnihilator([1,2,2;1,2,2],[1,2,2,4,5,6;3,2,1,4,5,6])
ans =
-2.8284 -2.8284 -2.1213 -5.6569 -7.0711 -8.4853
1.4142 0.0000 -0.7071 0.0000 0.0000 0.0000
>>
This routine zeroes out as many rows as possible in the rightmost block rows of the matrix(hmat).
"matlab/numericAnnihilateRows.m" 7 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function newhmat = numericAnnihilateRows(hmat)
[hrows,hcols]=size(hmat);
newhmat=numericComputeAnnihilator(hmat(:,hcols-hrows+1:hcols),hmat);
%newhmat=zapper*hmat;
<>
This routine zeroes out as many rows as possible in the leftmost block rows of the matrix(hmat).
"matlab/numericLeftAnnihilateRows.m" 8 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function newhmat = numericLeftAnnihilateRows(hmat)
[hrows,hcols]=size(hmat);
newhmat=numericComputeAnnihilator(hmat(:,1:hrows),hmat);
%newhmat=zapper*hmat;
<>
blue
AN EXAMPLE:
>> numericAnnihilateRows([2,4,3,9,2,2;7,3,4,1,2,2])
ans =
-6.3640 -4.9497 -4.9497 -7.0711 -2.8284 -2.8284
3.5355 -0.7071 0.7071 -5.6569 0.0000 0.0000
>>
>>
This routine implements Algorithm 1.
"matlab/numericRightAR.m" 9 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [auxConds,hstar]=numericRightAR(hmat)
fhmat=full(hmat);
oldAuxRows=-1;
[hrows,hcols]=size(hmat);
[auxConds,hstar]=numericShiftRightAndRecord(zeros(0,hcols-hrows),fhmat);
[newAuxRows,ignore]=size(auxConds);
while(oldAuxRows ~= newAuxRows)
oldAuxRows=newAuxRows;
[auxConds,hstar]=numericShiftRightAndRecord(auxConds,numericAnnihilateRows(hstar));
[newAuxRows,ignore]=size(auxConds);
end%end while
if(issparse(hmat))
auxConds=sparse(auxConds);
hstar=sparse(hstar);
end%if
1<>
blue
AN EXAMPLE:
>> hfv=sparse([0,0,-(1+0.1),0,1,1;0,-(1-0.6),0,1,0,0])
hfv =
(2,2) -0.4000
(1,3) -1.1000
(2,4) 1.0000
(1,5) 1.0000
(1,6) 1.0000
>>
>>
>> [zf,hf]=numericRightAR(hfv)
ans =
1
zf =
(1,2) -0.4000
(1,4) 1.0000
hf =
(1,3) 0.7778
(2,3) 0.7778
(1,4) 0.2828
(2,4) -0.2828
(1,5) -0.7071
(2,5) -0.7071
(1,6) -1.4142
(2,6) 0.0000
>>
>>
This routine implements Algorithm 1.
"matlab/numericLeftAR.m" 10 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [auxConds,hstar]=numericLeftAR(hmat)
fhmat=full(hmat);
oldAuxRows=-1;
[hrows,hcols]=size(hmat);
[auxConds,hstar]=numericShiftLeftAndRecord(zeros(0,hcols-hrows),fhmat);
[newAuxRows,ignore]=size(auxConds);
while(oldAuxRows ~= newAuxRows)
oldAuxRows=newAuxRows;
[auxConds,hstar]=numericShiftLeftAndRecord(auxConds,numericLeftAnnihilateRows(hstar));
[newAuxRows,ignore]=size(auxConds);
end%end while
if(issparse(hmat))
auxConds=sparse(auxConds);
hstar=sparse(hstar);
end%if
<>
blue
AN EXAMPLE:
hfv=sparse([0,0,-(1+0.1),0,1,1;0,-(1-0.6),0,1,0,0])
[zf,hf]=numericLeftAR(hfv)
hfv=sparse([0,0,-(1+0.1),0,1,1;0,-(1-0.6),0,1,0,0])
hfv =
(2,2) -0.4000
(1,3) -1.1000
(2,4) 1.0000
(1,5) 1.0000
(1,6) 1.0000
>>
>>
>> [zb,hb]=numericLeftAR(hfv)
[zb,hb]=numericLeftAR(hfv)
zb =
(1,1) -1.1000
(1,3) 1.0000
(1,4) 1.0000
hb =
(1,1) 1.1000
(2,2) -0.4000
(1,3) -1.0000
(1,4) -1.0000
(2,4) 1.0000
>>
>>
>>
>>
This routine applies Algorithm 1 in the forward and the reverse directions to provide inputs for the state space reduction routines. The auxiliary initial conditions span the invariant space associated with zero eigenvalues. Obtaining these vectors makes it possible to compute the minimal dimension transition matrix.
"matlab/numericBiDirectionalAR.m" 11 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [zb,hb,zf,hf]=numericBiDirectionalAR(hmat)
[hrows,hcols]=size(hmat);
[zf,hf]=numericRightAR(hmat);
'done forward'
%p=hcols:-1:1;
%ps=(hcols-hrows):-1:1;
%hrev=hmat(:,p);
[zb,hb]=numericLeftAR(hmat);
'done backward'
%[zbrows,zbcols]=size(zb);
%if zbrows > 0
%zb=zb(:,ps);
%end
%hb=hb(:,p);
<>
This routine implements Algorithm 2.
"matlab/numericAsymptoticConstraint.m" 12 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function qmat=numericAsymptoticConstraint(af,ab,hf,nlag)
[hrows,hcols]=size(hf);
[frows,ignore]=size(af);
[brows,ignore]=size(ab);
tm=numericTransitionMatrix(hf);
[j0,pimat,lila,qmat]=numericSqueeze([ab;af],hf);
[nrows,ncols]=size(lila);
[vm,lm]=numericSelectLargeIS(lila,min(nrows,(hcols-hrows)-frows-nlag*hrows));
bigOnes=numericEvExtend(j0,pimat,lm,vm,qmat);
%[tmlil,tmia]=numericEliminateInessentialLags(hf);
%[nrows,ncols]=size(tmlil);
%[vm,lm]=numericSelectLargeIS(tmlil,max(nrows,(hcols-hrows)-frows-nlag*hrows+2));
%[nrvm,ncvm]=size(vm);
%bigOnes(1:nrvm,tmia)=vm;
[nrb,nrc]=size(bigOnes);
qmat=[af;bigOnes];
<>
This routine implements Algorithm 2.
"matlab/numericAsymptoticConstraintContinuous.m" 13 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function qmat=numericAsymptoticConstraintContinuous(af,ab,hf,nlag)
[hrows,hcols]=size(hf);
[frows,ignore]=size(af);
[brows,ignore]=size(ab);
tm=numericTransitionMatrix(hf);
[j0,pimat,lila,qmat]=numericSqueeze([ab;af],hf);
[nrows,ncols]=size(lila);
[vm,lm]=numericSelectPositiveIS(lila,min(nrows,(hcols-hrows)-frows-nlag*hrows));
bigOnes=numericEvExtend(j0,pimat,lm,vm,qmat);
%[tmlil,tmia]=numericEliminateInessentialLags(hf);
%[nrows,ncols]=size(tmlil);
%[vm,lm]=numericSelectPositiveIS(tmlil,max(nrows,(hcols-hrows)-frows-nlag*hrows+2));
%[nrvm,ncvm]=size(vm);
%bigOnes(1:nrvm,tmia)=vm;
[nrb,nrc]=size(bigOnes);
qmat=[af;bigOnes];
<>
This routine implements Algorithm 3.
"matlab/numericCheckUniqueness.m" 14 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [bmat,infostring]=numericCheckUniqueness(hf,nlags,qmat)
[hrows,hcols]=size(hf);
[qrows,qcols]=size(qmat);
nleads=(qcols/hrows)-nlags;
if(qrows<nleads*hrows)
infostring='multiplicity of stable solutions'
bmat=[];
elseif(qrows>nleads*hrows)
infostring='no stable solutions except for particular initial conditions'
bmat=[];
else
bmat=numericAsymptoticAR(qmat);
infostring='unique stable solutions';
end
<>
blue
AN EXAMPLE:
[zb,hb,zf,hf]=numericBiDirectionalAR(hfv)
ans =
1
ans =
done forward
ans =
done backward
zb =
(1,1) -1.1000
(1,3) 1.0000
(1,4) 1.0000
hb =
(1,1) 1.1000
(2,2) -0.4000
(1,3) -1.0000
(1,4) -1.0000
(2,4) 1.0000
zf =
(1,2) -0.4000
(1,4) 1.0000
hf =
(1,3) 0.7778
(2,3) 0.7778
(1,4) 0.2828
(2,4) -0.2828
(1,5) -0.7071
(2,5) -0.7071
(1,6) -1.4142
(2,6) 0.0000
>> theq=numericAsymptoticConstraint(zf,zb,hf,1)
iter =
1
eigs =
1.1000
0.4000
stopcrit =
2.4728e-16
iter =
2
eigs =
1.1000
0.4000
stopcrit =
2.9673e-16
evrows =
1
evcols =
1
kronMat =
(1,1) 0.0818
(2,1) 0.9908
(1,2) 0.8343
(2,2) 0.0022
ans =
(1,1) 0.0818
(2,1) 0.9908
(1,2) 0.8343
(2,2) 0.0022
theq =
(2,1) 0.0000
(1,2) -0.4000
(2,2) -0.0000
(2,3) 1.2095
(1,4) 1.0000
(2,4) -0.6911
>>
>>
[zb,hb,zf,hf]=numericBiDirectionalAR(hfv)
theq=numericAsymptoticConstraint(zf,zb,hf,1)
LU factorization of (A-sigma*I) for sigma = 0.000000e+00: condest(U) = 3.285580e+00
iter =
1
eigs =
0.4000
1.1000
stopcrit =
2.4728e-16
iter =
2
eigs =
0.4000
1.1000
stopcrit =
2.4728e-16
theq =
(2,1) 0.0000
(1,2) -0.4000
(2,2) -0.0000
(2,3) 1.2095
(1,4) 1.0000
(2,4) -0.6911
>> [bnow,is]=numericCheckUniqueness(hfv,1,theq)
bnow =
(1,1) -0.0000
(1,2) 0.2286
(2,2) 0.4000
is =
unique stable solutions
>>
This routine uses the auxiliary initial conditions and the left invariant space vectors to construct an autoregressive representation of the stable solutions for the linear system.
"matlab/numericAsymptoticAR.m" 15 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function bmat=numericAsymptoticAR(qmat)
[qrows,qcols]=size(qmat);
bmat=-qmat(:,qcols-qrows+1:qcols)\qmat(:,1:qcols-qrows);
<>
The algorithm constructs transition matrices that are sparse and have many zero eigenvalues. The routines in this section compute a similarity transformation that allows subsequent routines to compute invariant space vectors in a space of lower dimension.
where
"matlab/numericEliminateInessentialLags.m" 16 =
function [a,js] = numericEliminateInessentialLags(h)
[neq,ncols]=size(h);
qcols=ncols-neq;
% [a,js] = build_a(h,qcols,neq)
%
% Build the companion matrix, deleting inessential lags.
% Solve for x_{t+nlead} in terms of x_{t+nlag},...,x_{t+nlead-1}.
left = 1:qcols;
right = qcols+1:qcols+neq;
hs=sparse(h);
%size(h(:,right))
%size(h(:,left))
hs(:,left) = -hs(:,right)\hs(:,left);
%h(:,left) = -h(:,right)\h(:,left);
% Build the big transition matrix.
a = zeros(qcols,qcols);
if(qcols > neq)
eyerows = 1:qcols-neq;
eyecols = neq+1:qcols;
a(eyerows,eyecols) = eye(qcols-neq);
end
hrows = qcols-neq+1:qcols;
a(hrows,:) = hs(:,left);
% Delete inessential lags and build index array js. js indexes the
% columns in the big transition matrix that correspond to the
% essential lags in the model. They are the columns of q that will
% get the unstable left eigenvectors.
%a=sparse(a);
js = 1:qcols;
zerocols = sum(abs(a)) == 0;
while( any(zerocols) )
a(:,zerocols) = [];
a(zerocols,:) = [];
js(zerocols) = [];
zerocols = sum(abs(a)) == 0;
end
ia = length(js);
a=sparse(a);
<>
This routine computes the 
components of the
mapping from reduced dimension eigenvectors to full dimensions eigenvectors.
"matlab/numericSqueeze.m" 17 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [j0,pi,a,qout]=numericSqueeze(auxFAuxB,preTransMat)
zeroTol=10^(-10);
%[tmlil,tmia]=numericEliminateInessentialLags(preTransMat);
tm=numericTransitionMatrix(preTransMat);
if(isempty(auxFAuxB))
afab=[];
notafab=speye(auxcols);
else
[auxRows,auxCols]=size(auxFAuxB);
[afab,notafab]=numericExtendToBasis(auxFAuxB);
end;%end if
qout=[afab;notafab];
if(isempty(auxFAuxB))
j0=[];
pi=[];
a=tm;
else
j0=afab * tm * afab.';
pi=notafab * tm * afab.';
a=notafab * tm * notafab.';
end;
<>
This routine computes the matrix for transforming reduced dimension eigenvectors into full dimension eigenvectors.
"matlab/numericMapper.m" 18 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function kronMat=numericMapper(bigJ0,bigPi,evMat)
[jrows,jcols]=size(bigJ0);
[evrows,evcols]=size(evMat)
kronMat=(kron(speye(jrows),evMat) - kron(bigJ0.',speye(evrows)))\kron(bigPi.',speye(evrows))
<>
"matlab/numericEvExtend.m" 19 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function bigEvecs=numericEvExtend(bigJ0,bigPi,evMat,yMat,qmat)
[jrows,jcols]=size(bigJ0);
[evrows,evcols]=size(evMat);
[yrows,ycols]=size(yMat);
numericMapper(bigJ0,bigPi,evMat)
if(jrows>0)
kronMat=(kron(speye(jrows),evMat) - kron(bigJ0.',speye(evrows)))\kron(bigPi.',speye(evrows));
xMat=reshape(kronMat * reshape(yMat,yrows*ycols,1),yrows,jcols);
bigEvecs=[xMat,yMat]*qmat;
else
bigEvecs=yMat*qmat;
end;
<>
This routine identifies the eigenvalues with magnitude larger than one. The routine returns a matrix whose rows span the left invariant space of the transition matrix, along with the corresponding diagonal matrix of eigenvalues.
"matlab/numericSelectLargeIS.m" 20 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [bigv,biglam]=numericSelectLargeIS(tmat,numToComp)
[allv,alllam]=eigs(tmat.',numToComp);
[iv,jv]=find(abs(alllam)>1);
bigv=allv(:,iv).';
biglam=alllam(iv,jv);
<>
This routine identifies the positive eigenvalues. The routine returns a matrix whose rows span the left invariant space of the transition matrix, along with the corresponding diagonal matrix of eigenvalues.
"matlab/numericSelectPositiveIS.m" 21 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [bigv,biglam]=numericSelectPositiveIS(tmat,numToComp)
[allv,alllam]=eigs(tmat.',numToComp,'LR');
[iv,jv]=find(alllam>0);
bigv=allv(:,iv).';
biglam=alllam(iv,jv);
<>
"matlab/numericExtendToBasis.m" 22 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [orig,extension]=numericExtendToBasis(matRows)
zeroTol=10^(-10);
[nrows,ncols]=size(matRows);
v=orth([full(matRows);eye(ncols)]');
orig=v(:,1:nrows)';%assuming full rank
extension=(v(:,nrows+1:ncols))';%assuming full rank
<>
Following [[]],
For time t expectations, one has
One can construct a linear transition matrix for iterating the system forward
One can exploit the relation 
to compute
the unconditional variance( 
):
This routine computes the time t observable structure.
The routine returns the routine returns the S matrix and the
component of 
that multiplies the 
, 
.
"matlab/numericObservableStructure.m" 23 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [smat,sjac]=numericObservableStructure(hmat,qmat)
[hrows,hcols]=size(hmat);
[qrows,qcols]=size(qmat);
smat=[hmat,-eye(hrows);[zeros(qrows,hcols-qcols),qmat,zeros(qrows,hrows)]];
sjac=inv(smat(:,hcols+hrows-(hrows+qrows)+1:hcols+hrows));
smat=sjac *smat;
smat=smat(qrows+1:qrows+hrows,1:hcols-qrows);
sjac=smat(:,hcols-qrows-hrows+1:hcols-qrows);
<>
This routine computes the time t-1 observable structure.
The routine returns the routine returns the S matrix and the
component of 
that multiplies the 
, 
.
For t-1 expectations we have
"matlab/numericObservableStructureTM1.m" 24 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [smat,sjac]=numericObservableStructureTM1(hmat,bmat)
[hrows,hcols]=size(hmat);
[brows,bcols]=size(bmat);
smat= hmat(:,bcols+1:bcols+hrows) * [-bmat,eye(hrows)];
sjac=inv(smat(:,bcols+1:bcols+hrows));
<>
blue
AN EXAMPLE:
blue
AN EXAMPLE:
[thes,thesinv]=numericObservableStructure(hfv,theq)
[thes,thesinv]=numericObservableStructure(hfv,theq)
thes =
(2,2) 0.4000
(1,3) 1.1000
(1,4) -0.6286
(2,4) -1.0000
thesinv =
(1,1) 1.1000
(1,2) -0.6286
(2,2) -1.0000
>>
>>
[thesTM1,thesinvTM1]=numericObservableStructureTM1(hfv,bnow)
thesTM1 =
(1,1) -0.0000
(1,2) 0.2514
(2,2) -0.4000
(1,3) -1.1000
(2,4) 1.0000
thesinvTM1 =
(1,1) -0.9091
(2,2) 1.0000
>>
"matlab/numericUnconditionalCovariance.m" 25 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function varCov=numericUnconditionalCovariance(smat,shockCov)
[srows,scols]=size(smat);
sinv=inv(smat(:,scols-srows+1:scols));
varTrans=[zeros(scols-(2*srows),srows),speye(scols-(2*srows));-sinv*smat(:,1:scols-srows)];
varCov=(speye((scols-srows)^2)-kron(varTrans,varTrans)) \ reshape([zeros(scols-(2*srows),scols-srows);zeros(srows,scols-(2*srows)),sinv*shockCov*sinv'],(scols-srows)^2,1);
<>
This routine applies the above routines to transform a linear model to produce an unconstrained autoregressive representation a constrained autoregression and observable structure
"matlab/numericAim.m" 26 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
function [q,b,info,sinv]=numericAim(nlag,hmat)
[ab,hb,af,hf]=numericBiDirectionalAR(hmat);
'done biDir'
q=numericAsymptoticConstraint(af,ab,hf,nlag);
'done q comp'
[b,info]=numericCheckUniqueness(hf,nlag,q);
'done b comp'
[s,sinv]=numericObservableStructure(hmat,q);
<>
blue
AN EXAMPLE:
[qm,bm,im,sm]=numericAim(1,hfv)
ans =
1
ans =
done forward
ans =
done backward
ans =
done biDir
iter =
1
eigs =
1.1000
0.4000
stopcrit =
1.9782e-16
iter =
2
eigs =
1.1000
0.4000
stopcrit =
1.4837e-16
evrows =
1
evcols =
1
kronMat =
(1,1) 0.0818
(2,1) 0.9908
(1,2) 0.8343
(2,2) 0.0022
ans =
(1,1) 0.0818
(2,1) 0.9908
(1,2) 0.8343
(2,2) 0.0022
ans =
done q comp
ans =
done b comp
qm =
(2,1) 0.0000
(1,2) -0.4000
(2,2) -0.0000
(2,3) 1.2095
(1,4) 1.0000
(2,4) -0.6911
bm =
(1,1) -0.0000
(1,2) 0.2286
(2,2) 0.4000
im =
unique stable solutions
sm =
(1,1) 1.1000
(1,2) -0.6286
(2,2) -1.0000
>>
"matlab/full.m" 27 =
function fl=full(mat)
fl=mat;
<>
"matlab/sparse.m" 28 =
function fl=sparse(mat)
fl=mat;
<>
"matlab/issparse.m" 29 =
function res=issparse(mat)
res=0;
<>
"matlab/numericAIMVersion.m" 30 =
function res=numericAIMVersion()
res='$Revision: 1.14 $ $Date: 1998/09/09 18:34:43 $'
<>
"numericLinearAimTestMLab.m" 31 =
%$Author: m1gsa00 $
%$Date: 1998/09/09 18:34:43 $
%$Id: reliableMLab.w,v 1.14 1998/09/09 18:34:43 m1gsa00 Exp m1gsa00 $
path(path,'/irm/home/m1gsa00/aim/reliable/matlab')
hfv=sparse([0,0,-(1+0.1),0,1,1;0,-(1-0.6),0,1,0,0])
jefftest=[-0.9000, 0, 1.0000, 0, 0, 0; 0 0 , 0 , 1.0000 , -0.9000 , 0]
<>
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
blue
AN EXAMPLE:
blue
AN EXAMPLE:
blue
AN EXAMPLE:
blue
AN EXAMPLE:
blue
AN EXAMPLE:
blue
AN EXAMPLE:
blue
AN EXAMPLE:
blue
AN EXAMPLE:
hfv =
(2,2) -0.4000
(1,3) -1.1000
(2,4) 1.0000
(1,5) 1.0000
(1,6) 1.0000
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 32 =
[ab,hb,af,hf]=numericBiDirectionalAR(hfv)
[atb,htb,atf,htf]=numericBiDirectionalAR(jefftest)
<>
ab =
(1,1) -1.1000
(1,3) 1.0000
(1,4) 1.0000
hb =
(2,1) -1.1000
(1,2) -0.4000
(2,3) 1.0000
(1,4) 1.0000
(2,4) 1.0000
af =
(1,2) -0.4000
(1,4) 1.0000
hf =
(1,3) -1.1000
(2,4) -0.4000
(1,5) 1.0000
(1,6) 1.0000
(2,6) 1.0000
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 33 =
tm=numericTransitionMatrix(hf)
<>
tm =
(1,3) 1.0000
(3,3) 1.1000
(2,4) 1.0000
(3,4) -0.4000
(4,4) 0.4000
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 34 =
[j0,pimat,lila,qout]=numericSqueeze([ab;af],hf)
[jt0,pimatt,lilat,qlout]=numericSqueeze([atb;atf],htf)
<>
j0 =
1.0e-16 *
(1,1) -0.3098
(1,2) -0.4149
pimat =
(1,1) -0.8017
(2,1) -0.4275
(1,2) 0.4858
(2,2) -0.9882
lila =
(1,1) 0.9667
(2,1) -0.1511
(1,2) -0.5001
(2,2) 0.5333
qu =
(1,1) 0.6140
(2,1) -0.3720
(2,2) 0.4343
(1,3) -0.5581
(2,3) 0.3382
(1,4) -0.5581
(2,4) -0.7474
ql =
(1,1) 0.6962
(1,2) 0.2321
(2,2) 0.8704
(1,3) 0.6730
(2,3) -0.3482
(1,4) 0.0928
(2,4) 0.3482
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 35 =
[vbig,lambig]=eigs(tm.');
[vlil,lamlil]=eigs(lila.');
%[nea,neb]=numericExtendToBasis([ab;af]);
bigOnes=numericEvExtend(j0,pimat,lamlil,vlil,qout)
<>
bigOnes =
-0.0000 0 1.2095 -0.6911
0.0000 -0.0000 -0.0000 2.7753
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 36 =
[vm,lm]=numericSelectLargeIS(tm,2)
<>
vm =
0.0000 0.0000 -0.8682 0.4961
lm =
1.1000
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 37 =
theq=numericAsymptoticConstraint(af,ab,hf,1)
<>
theq =
(1,2) -0.4000
(2,2) 0.0000
(2,3) 1.2095
(1,4) 1.0000
(2,4) -0.6911
File defined by scraps 31, 32, 33, 34, 35, 36, 37, 38.
"numericLinearAimTestMLab.m" 38 =
theb=numericAsymptoticAR(theq)
[thes,thesinv]=numericObservableStructure(hfv,theq)
covm= numericUnconditionalCovariance(thes,speye(2))
<>
theb =
(1,1) -0.0000
(1,2) -0.2286
(2,2) -0.4000
[thes,thesinv]=numericObservableStructure(hfv,theq)
covm= numericUnconditionalCovariance(thes,speye(2))
[thes,thesinv]=numericObservableStructure(hfv,theq)
thes =
(2,2) 0.4000
(1,3) 1.1000
(1,4) -0.6286
(2,4) -1.0000
thesinv =
(1,1) 1.1000
(1,2) -0.6286
(2,2) -1.0000
>>
>> covm= numericUnconditionalCovariance(thes,speye(2))
covm =
(1,1) 1.2152
(2,1) 0.6803
(3,1) 0.6803
(4,1) 1.1905
>>
>>
"matlab/full.m"
"matlab/issparse.m"
"matlab/numericAim.m"
"matlab/numericAIMVersion.m"
"matlab/numericAnnihilateRows.m"
"matlab/numericAsymptoticAR.m"
"matlab/numericAsymptoticConstraint.m"
"matlab/numericAsymptoticConstraintContinuous.m"
"matlab/numericBiDirectionalAR.m"
"matlab/numericCheckUniqueness.m"
"matlab/numericComputeAnnihilator.m"
"matlab/numericEliminateInessentialLags.m"
"matlab/numericEvExtend.m"
"matlab/numericExtendToBasis.m"
"matlab/numericLeftAnnihilateRows.m"
"matlab/numericLeftAR.m"
"matlab/numericLeftMostAllZeroQ.m"
"matlab/numericMapper.m"
"matlab/numericObservableStructure.m"
"matlab/numericObservableStructureTM1.m"
"matlab/numericRightAR.m"
"matlab/numericRightMostAllZeroQ.m"
"matlab/numericSelectLargeIS.m"
"matlab/numericSelectPositiveIS.m"
"matlab/numericShiftLeftAndRecord.m"
"matlab/numericShiftRightAndRecord.m"
"matlab/numericSqueeze.m"
"matlab/numericTransitionMatrix.m"
"matlab/numericUnconditionalCovariance.m"
"matlab/sparse.m"
"numericLinearAimTestMLab.m"
numericAnnihilateRows:
numericBiDirectionalAR:
numericComputeAnnihilator:
numericEvExtend:
numericExtendToBasis:
numericLeftAR:
numericLeftMostAllZeroQ:
numericMapper:
numericRightAR:
numericRightMostAllZeroQ:
numericSelectLargeIS:
numericSelectPositiveIS:
numericShiftLeftAndRecord:
numericShiftRightAndRecord:
numericSqueeze:
numericTransitionMatrix:
zeroTol:
.
[Anderson]Anderson1997bANDER:PARA Anderson, Gary (1997b). A parallel programming implementation of the anderson-moore(aim) algorithm. Unpublished Manuscript, Board of Governors of the Federal Reserve System. Downloadable copies of this and other related papers at http://www.federalreserve.gov/pubs/oss/oss4/aimindex.html.
[Anderson and Moore]Anderson and Moore1983ANDER:AIM1 Anderson, Gary and George Moore (1983). An efficient procedure for solving linear perfect foresight models. Unpublished Manuscript, Board of Governors of the Federal Reserve System. Downloadable copies of this and other related papers at http://www.federalreserve.gov/pubs/oss/oss4/aimindex.html.
[Anderson and Moore]Anderson and Moore1985ANDER:AIM2 Anderson, Gary and George Moore (1985). A linear algebraic procedure for solving linear perfect foresight models. Economics Letters.
[Berry]Berry1996berry Berry, Michael W. (1996). Large scale sparse singular value computations. University of Tennessee, Department of Computer Science.
[Blanchard and Kahn]Blanchard and Kahn1980blanchard80 Blanchard, Olivier Jean and C. Kahn (1980). The solution of linear difference models under rational expectations. Econometrica.
[Fuhrer]Fuhrer1994fuhrer94 Fuhrer, Jeff (1994). Optimal monetary policy in a model of overlapping price contracts. Technical Report 94-2. Federal Reserve Bank of Boston.
[Golub and van Loan]Golub and van Loan1989golub89 Golub, Gene H. and Charles F. van Loan (1989). Matrix Computations. Johns Hopkins.
[Krishnamurthy]Krishnamurthy1989krishnamurthy89 Krishnamurthy, E. V. (1989). Parallel Processing: Principles and Practice. Addison-Wesley.
[Noble]Noble1969NOBLE Noble, Ben (1969). Applied Linear Algebra. Prentice-Hall, Inc.
[Taylor]Taylor77taylor77 Taylor, J. (77). Conditions for unique solutions in stochastic macroeconomic models with rational expectations. Econometrica 45, 1377-1385.
[Whiteman]Whiteman1983whiteman83 Whiteman, C. H. (1983). Linear Rational Expectations Models: A User's Guide. University of Minnesota.