AMA to AndersonMoore

src/AndersonMoore.jl

@@ -0,0 +1,21 @@
+__precompile__()
+
+module AndersonMoore
+# http://www.stochasticlifestyle.com/finalizing-julia-package-documentation-testing-coverage-publishing/
+
+# Set-up for callSparseAim
+const lib = Libdl.dlopen(normpath(joinpath(dirname(@__FILE__), "..", "deps", "libAndersonMoore")))
+const sym = Libdl.dlsym(lib, :callSparseAim)
+
+# Include all files    
+for (root, dirs, files) in walkdir(joinpath(dirname(@__FILE__)))
+    for file in files
+        file == "AndersonMoore.jl" ? nothing : include(joinpath(root, file))
+    end
+end
+
+# Export all functions
+export exactShift!, numericShift!, shiftRight!, buildA!, augmentQ!, eigenSys!, reducedForm,
+AndersonMooreAlg, sameSpan, deleteCols, deleteRows, callSparseAim, checkAM, err, gensysToAMA
+
+end # module

src/AndersonMoore/AndersonMooreAlg.jl

@@ -0,0 +1,108 @@
+"""
+    AndersonMooreAlg(h, qcols, neq)
+
+Solve a linear perfect foresight model using the julia eig
+function to find the invariant subspace associated with the big
+roots.  This procedure will fail if the companion matrix is
+defective and does not have a linearly independent set of
+eigenvectors associated with the big roots.
+
+  Input arguments:
+
+    h         Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
+    neq       Number of equations.
+    nlag      Number of lags.
+    nlead     Number of leads.
+    condn     Zero tolerance used as a condition number test
+              by numericShift and reducedForm.
+    upper     Inclusive upper bound for the modulus of roots
+              allowed in the reduced form.
+ 
+  Output arguments:
+ 
+    b         Reduced form coefficient matrix (neq,neq*nlag).
+    rts       Roots returned by eig.
+    ia        Dimension of companion matrix (number of non-trivial
+              elements in rts).
+    nexact    Number of exact shiftRights.
+    nnumeric  Number of numeric shiftRights.
+    lgroots   Number of roots greater in modulus than upper.
+    AMAcode   Return code: see function AMAerr.
+"""
+function AndersonMooreAlg(hh::Array{Float64,2}, neq::Int64, nlag::Int64, nlead::Int64, anEpsi::Float64, upper::Float64) 
+
+    if(nlag < 1 || nlead < 1) 
+        error("AMA_eig: model must have at least one lag and one lead.")
+    end
+
+    # Initialization.
+    nexact   = 0
+    nnumeric = 0
+    lgroots  = 0
+    iq       = 0
+    AMAcode  = 0
+    bb       = 0
+    qrows    = neq * nlead
+    qcols    = neq * (nlag + nlead)
+    bcols    = neq * nlag
+    qq       = zeros(qrows, qcols)
+    rts      = zeros(qcols, 1)
+
+    # Compute the auxiliary initial conditions and store them in q.
+
+    (hh, qq, iq, nexact) = exactShift!(hh, qq, iq, qrows, qcols, neq)
+    if (iq > qrows) 
+        AMAcode = 61
+        return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
+    end
+
+    (hh, qq, iq, nnumeric) = numericShift!(hh, qq, iq, qrows, qcols, neq, anEpsi)
+    if (iq > qrows) 
+        AMAcode = 62
+        return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
+    end
+
+    #  Build the companion matrix.  Compute the stability conditions, and
+    #  combine them with the auxiliary initial conditions in q.  
+
+    (aa, ia, js) = buildA!(hh, qcols, neq)
+
+    if (ia != 0)
+        for element in aa
+            if isnan(element) || isinf(element)
+                display("A is NAN or INF")
+                AMAcode = 63
+                return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
+            end
+        end
+        (ww, rts, lgroots) = eigenSys!(aa, upper, min(size(js, 1), qrows - iq + 1))
+
+
+        qq = augmentQ!(qq, ww, js, iq, qrows)
+    end
+
+    test = nexact + nnumeric + lgroots
+    if (test > qrows)
+        AMAcode = 3
+    elseif (test < qrows)
+        AMAcode = 4
+    end
+
+    # If the right-hand block of q is invertible, compute the reduced form.
+
+    if(AMAcode == 0)
+        (nonsing,bb) = reducedForm(qq, qrows, qcols, bcols, neq, anEpsi)
+        if ( nonsing && AMAcode==0)
+            AMAcode =  1
+        elseif (!nonsing && AMAcode==0)
+            AMAcode =  5
+        elseif (!nonsing && AMAcode==3)
+            AMAcode = 35
+        elseif (!nonsing && AMAcode==4)
+            AMAcode = 45
+        end
+    end
+
+    return bb, rts, ia, nexact, nnumeric, lgroots, AMAcode
+ # (bbJulia,rtsJulia,iaJulia,nexJulia,nnumJulia,lgrtsJulia,AMAcodeJulia) = AMAalg(hh,neq,nlag,nlead,anEpsi,1+anEpsi)
+end # AndersonMooreAlg

src/AndersonMoore/augmentQ!.jl

@@ -0,0 +1,18 @@
+"""
+    augmentQ(qq, ww, js, iq, qrows)
+
+Copy the eigenvectors corresponding to the largest roots into the
+remaining empty rows and columns js of q 
+"""
+function augmentQ!(qq::Array{Float64,2}, ww::Array{Float64,2}, js::Array{Int64,2}, iq::Int64, qrows::Int64) 
+
+    if(iq < qrows)
+        lastrows = (iq + 1) : qrows
+        wrows    = 1 : length(lastrows)
+        qq[lastrows, js] = @views ww[:, wrows]'
+    end
+
+    return qq
+
+end # augmentQ
+

src/AndersonMoore/buildA!.jl

@@ -0,0 +1,57 @@
+"""
+    buildA(hh, qcols, neq)
+
+Build the companion matrix, deleting inessential lags.
+Solve for x_{t+nlead} in terms of x_{t+nlag},...,x_{t+nlead-1}.
+
+"""
+function buildA!(hh::Array{Float64,2}, qcols::Int64, neq::Int64) 
+
+    left  = 1:qcols
+    right = (qcols + 1):(qcols + neq)
+
+    tmp = similar(hh)
+    tmp[:, left] =  \(-hh[:, right], hh[:, left])
+
+    #  Build the big transition matrix.
+
+    aa = zeros(qcols, qcols)
+
+    if(qcols > neq)
+        eyerows = 1:(qcols - neq)
+        eyecols = (neq + 1):qcols
+        aa[eyerows, eyecols] = eye(qcols - neq)
+    end
+    hrows      = (qcols - neq + 1):qcols
+    aa[hrows, :] =  tmp[:, 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. 
+
+    js       = Array{Int64}(1,qcols)
+    for ii in 1 : qcols
+        js[1, ii] = ii
+    end
+
+    zerocols = sum(abs.(aa), 1)
+    zerocols = find(col->(col == 0), zerocols)
+
+
+     while length(zerocols) != 0
+        # aa = filter!(x->(x !in zerocols), aa)
+
+        aa = deleteCols(aa, zerocols)        
+        aa = deleteRows(aa, zerocols)
+        js = deleteCols(js, zerocols)
+
+        zerocols = sum(abs.(aa), 1)  
+        zerocols = find(col->(col == 0), zerocols)
+
+    end
+    ia = length(js)
+
+    return (aa::Array{Float64,2}, ia::Int64, Int64.(js)::Array{Int64,2})
+
+end # buildA

src/AndersonMoore/eigenSys!.jl

@@ -0,0 +1,34 @@
+"""
+    eigenSys!(h, qcols, neq)
+
+Compute the roots and the left eigenvectors of the companion
+matrix, sort the roots from large-to-small, and sort the
+eigenvectors conformably.  Map the eigenvectors into the real
+domain. Count the roots bigger than uprbnd.
+"""
+function eigenSys!(aa::Array{Float64,2}, upperbound::Float64, rowsLeft::Int64) 
+
+    roots = eigvals(aa')
+    ww = eigvecs(aa')
+
+    # sort eigenvalues in descending order of magnitude
+    magnitudes = abs.(roots)
+    highestToLowestMag = sortperm(magnitudes, rev = true) # reverse order
+    roots = roots[highestToLowestMag]
+
+    # sort eigenvecs in order found above
+    ww = ww[:, highestToLowestMag]
+
+    #  Given a complex conjugate pair of vectors W = [w1,w2], there is a
+    #  nonsingular matrix D such that W*D = real(W) + imag(W).  That is to
+    #  say, W and real(W)+imag(W) span the same subspace, which is all
+    #  that AMA cares about. 
+
+    ww = ( real(ww) + imag(ww) )::Array{Float64,2}
+
+    # count how many roots are above the upperbound threshold
+    lgroots = mapreduce((mag->mag > upperbound ? 1 : 0), +, magnitudes)
+
+    return (ww, roots, lgroots)
+
+end # eigenSys!

src/AndersonMoore/err.jl

@@ -0,0 +1,46 @@
+function err(c) 
+# err(c)
+#
+# Interpret the return codes generated by the AMA routines.
+#
+# The return code c = 2 is used by AMA_schur.m but not by AMA_eig.m.
+
+# Original author: Gary Anderson
+# Original file downloaded from:
+# http://www.federalreserve.gov/Pubs/oss/oss4/code.html
+# Adapted for Dynare by Dynare Team.
+#
+# This code is in the public domain and may be used freely.
+# However the authors would appreciate acknowledgement of the source by
+# citation of any of the following papers:
+#
+# Anderson, G. and Moore, G.
+# "A Linear Algebraic Procedure for Solving Linear Perfect Foresight
+# Models."
+# Economics Letters, 17, 1985.
+#
+# Anderson, G.
+# "Solving Linear Rational Expectations Models: A Horse Race"
+# Computational Economics, 2008, vol. 31, issue 2, pages 95-113
+#
+# Anderson, G.
+# "A Reliable and Computationally Efficient Algorithm for Imposing the
+# Saddle Point Property in Dynamic Models"
+# Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
+# pages 472-489
+
+    if(c==1)  e="AndersonMoore: unique solution."
+elseif(c==2)  e="AndersonMoore: roots not correctly computed by real_schur."
+elseif(c==3)  e="AndersonMoore: too many big roots."
+elseif(c==35) e="AndersonMoore: too many big roots, and q(:,right) is singular."
+elseif(c==4)  e="AndersonMoore: too few big roots."
+elseif(c==45) e="AndersonMoore: too few big roots, and q(:,right) is singular."
+elseif(c==5)  e="AndersonMoore: q(:,right) is singular."
+elseif(c==61) e="AndersonMoore: too many exact shiftrights."
+elseif(c==62) e="AndersonMoore: too many numeric shiftrights."
+else          e="AndersonMoore: return code not properly specified"
+end
+
+return e
+
+end

src/AndersonMoore/exactShift!.jl

@@ -0,0 +1,52 @@
+"""
+    exactShift!(hh,qq, iq, qrows, qcols, neq)
+
+Compute the exact shiftrights and store them in qq.
+"""
+function exactShift!(hh::Array{Float64,2}, qq::Array{Float64,2}, iq::Int64, qRows::Int64, qCols::Int64, neq::Int64) 
+
+    # total number of shifts
+    nexact = 0
+
+    # functions to seperate hh
+    left = 1:qCols
+    right = (qCols + 1):(qCols + neq)
+
+
+    # get right most columns of hh
+    zerorows = hh[:, right]'    
+
+    # compute absolute value  
+    zerorows = abs.(zerorows)
+
+    # take the sum of the rows (aka 1st dimenison in julia)
+    zerorows = sum(zerorows, 1)
+
+    # anon function returns index of the rows who sum is 0
+    zerorows = find(row->(row == 0), zerorows)
+
+    # continues until all zerorow indexes have been processed
+    while length(zerorows) != 0 && (iq <= qRows)
+        nz = size(zerorows, 1)
+
+        # insert only the indexes found with zerorows into qq
+        qq[(iq + 1):(iq + nz), :] = hh[zerorows, left]
+
+	# update by shifting right by $neq columns
+        hh[zerorows,:] = shiftRight!(hh[zerorows, :], neq)
+
+        # increment our variables the amount of zerorows found
+        iq = iq + nz
+        nexact = nexact + nz
+
+        # update zerorows as before but with our new hh matrix
+        zerorows = hh[:, right]'
+        zerorows = abs.(zerorows)
+        zerorows = sum(zerorows, 1)
+        zerorows = find(row->(row == 0), zerorows)
+    end # while
+
+    return (hh, qq, iq, nexact)  
+
+end # exactShift
+

src/AndersonMoore/numericShift!.jl

@@ -0,0 +1,48 @@
+"""
+    numericShift!(hh, qq, iq, qrows, qcols, neq, condn)
+
+Compute the numeric shiftrights and store them in q.
+"""
+function numericShift!(hh::Array{Float64,2}, qq::Array{Float64,2}, iq::Int64, qRows::Int64, qCols::Int64, neq::Int64, condn::Float64) 
+
+    # total number of shifts
+    nnumeric = 0
+
+    # functions to seperate hh
+    left = 1:qCols
+    right = (qCols + 1):(qCols + neq)
+
+    # preform QR factorization on right side of hh
+    F = qrfact(hh[:, right], Val{true})
+
+    # filter R only keeping rows that are zero
+    zerorows = abs.(diag(F[:R]))
+    zerorows = find(x->(float(x) <= condn), zerorows)
+
+    @inbounds while (length(zerorows) != 0) && (iq <= qRows)
+        # update hh with matrix multiplication of Q and hh
+        hh = *(F[:Q]', hh)
+
+        # need total number of zero rows
+        nz = size(zerorows, 1)
+
+        # update qq to keep track of rows that are zero
+        qq[(iq + 1):(iq + nz), :] = hh[zerorows, left]
+
+        # update hh by shifting right
+        hh[zerorows, :] = shiftRight!(hh[zerorows, :], neq)
+
+        # increment our variables by number of shifts
+        iq = iq + nz
+        nnumeric = nnumeric + nz
+
+        # redo QR factorization and filter R as before
+        F = qrfact(hh[:, right], Val{true})
+        zerorows = abs.(diag(F[:R]))
+        zerorows = find(x->(float(x) <= condn), zerorows)
+
+    end # while
+
+    return(hh, qq, iq, nnumeric)
+
+end # numericShift