[Python] Givens Rotation on Python

I’m replicating Shapiro & Watson (1988) and need to impose the long-run restriction with a Givens rotation. The target identity is C_inf@B=C_inf(P@Q_total)=(C_inf@P)@Q_total=A@Q_total

Yet my code gives C(1)@B differnt from A@Q_total

I suspect I’ve defined or applied the Givens rotation incorrectly. Code is below—any guidance is welcome. def givens_rotation(a, b): r = np.hypot(a, b) if r == 0: return 1.0, 0.0 c = a / r s = b / r return c, s

def apply_givens(M, i, j, k):

x = M[i, j]
y = M[i, k]
c, s = givens_rotation(x, y)

n = M.shape[0]
G = np.eye(n)
G[j, j] = c
G[j, k] = s
G[k, j] = -s
G[k, k] = c

M[:] = M @ G
return G


zero_constraints = [
(0, 4), # step 1: zero A[0,4]
(0, 3), # step 2: zero A[0,3]
(0, 2), # step 3: zero A[0,2]
(0, 1), # step 4: zero A[0,1]
#
(1, 4), # step 5: zero A[1,4]
(1, 3), # step 6: zero A[1,3]
(1, 2), # step 7: zero A[1,2]
#
(2, 4), # step 8: zero A[2,4]
(2, 3) # step 9: zero A[2,3]
]


for (i, j) in zero_constraints:
k = n - 1
k = 4 if j < 4 else (3 if j < 3 else (2 if j < 2 else (1 if j < 1 else None)))

G_jk = apply_givens(A, i, j, k)

Continue reading...