Snippets

Matrix

• Create
• Kalman filter
• Eigenvalues

Create

Create a matrix

# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Create matrix

# Shows how to create a matrix
# and do simple matrix
# calculations

# A has 2 rows and 3 columns
A = matrix( 2, 3,   1, 2, 3,
                    4, 5, 6 );

# B has 3 rows and 3 columns
B = matrix( 3, 3,   1,  1,  2,
                    3,  5,  8,
                   13, 21, 34 );

# Calculate product C = A * B
# and inverse D = 1 / B
C = A * B;
D = 1 / B;

# Print output, see log for
# details
print( A, B, C, D );

Kalman filter

Calculate HP filter using KF

# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Kalman filter

# This calculates the HP
# filter using the Kalman
# filter
# The approximation is not
# exact but is good, especially
# for longer timer series

# Download US GDP from FRED
GDP = dataFred( "GDPC1" );

# We will apply the filter to
# log of GDP
lGDP = ln( GDP );

# HP filter parameter
lambda = 1600;
li     = 1 / lambda;

# Some auxilliary variables
j = firstDate( lGDP );
n = count    ( lGDP );
s = lGDP[ j ];

# Prepare the Kalman filter

# State transition F 2x2
F = matrix( 2, 2, 2, -1, 1, 0 );
# State variance Q 2x2
Q = matrix( 2, 2, li );

# Observation map H 2x1 and H'
H = matrix( 2, 1, 1, 0 );
Ht = matrixT( H );
# Observation variance
R = matrix( 1, 1, 1 );

# Prediction state
# matrixX repeats s throughout
Z0 = matrixX( s, 2, n );
P0 = matrix( 2, 2*n, 10000 );
# Prediction update
Z1 = matrix( 2,   n );
P1 = matrix( 2, 2*n );
# Backward pass
ZB = matrix( 2,   n );
PB = matrix( 2, 2*n );

# Simulate the Kalman filter
# (Just 1 iteration)
iter = 1;
while( iter <= 1,

  # Forward pass
  # Load initial state
  ZU = Z0[ 1:2; 1 ],
  PU = P0[ 1:2; 1:2 ],

  i = 1,
  while( i <= n,
    # Predict state
    ZP = F * ZU,
    PP = F * PU *
         matrixT( F ) + Q,

    # Kalman gain
    K = PP * H /
        ( Ht * PP * H + R ),

    # Update
    ZU = ZP + K *
         ( lGDP[ j + i - 1 ] -
           Ht * ZP ),
    PU = PP - K * Ht * PP,

    # Store
    Z0[ 1:2; i         ] = ZP,
    P0[ 1:2; i*2-1:i*2 ] = PP,
    Z1[ 1:2; i         ] = ZU,
    P1[ 1:2; i*2-1:i*2 ] = PU,

    i = i + 1 ),

  # Backward pass

  # Start with last update
  ZL = ZU,
  PL = PU,

  # Store
  ZB[ 1:2; n ] = ZL,
  PB[ 1:2; 2*n-1:2*n ] = PL,

  # Now execute backward pass
  i = n - 1,
  while( i >= 1,
    ZP = Z0[ 1:2; i+1 ],
    PP = P0[ 1:2; (i+1)*2-1:
                  (i+1)*2],
    ZU = Z1[ 1:2; i ],
    PU = P1[ 1:2; i*2-1:i*2],

    L = PU * matrixT( F ) / PP,

    ZL = ZU + L * ( ZL - ZP ),
    PL = PU + L * ( PL - PP ) *
         matrixT( L ),

    # Store
    ZB[ 1:2; i ] = ZL,
    PB[ 1:2; 2*i-1:2*i ] = PL,

    i = i - 1 ),

  iter = iter + 1 );

# Now HP filter is in ZB
print( matrixT( ZB ) );

# We call it 'potential'
lPot = ZB[ 1; 1:n ];
print( matrixT( lPot ) );
Pot = exp( lPot );

# Create an observation matrix
# with these values concatenated
obsMat = matrixCatColumns(
            matrixObs( GDP ),
            matrixT( Pot ) );

# Calculate gap
# Column 2 divide by column 3 - 1
gap = matrixAsSeries(
        obsMat[ 1:n; 1,2 ] ) /
      matrixAsSeries(
        obsMat[ 1:n; 1,3 ] ) - 1;
print( obsMat );
print( gap );

# Also calculate using
# considerably easier
# built in Kalman filter
# function 'kalmanFilter'
print( kalmanFilter( "-",
       matrix( 2, 2, 2, -1, 1, 0 ),
       matrix( 2, 2, 1/1600 ),
       lGDP, matrix( 2, 1, 1 ),
       matrix( 2, 1, 0 ),
       matrix( 2, 2, 10000 ) ) );

Eigenvalues

Calculate eigenvalues and eigenvectors

# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Eigenvalues

# Create a random matrix
A = matrixRand( 10, 10 ) - 0.5;

# Some formatting
format( 10 );
formatFix();

# Calculate the eigenvalues
# This returns two columns,
# one for the real and one
# for the imaginary part
L = matrixEig( A );

# Get the first eigenvalue and
# split it in its real and
# imaginary parts
L1 = matrixGetRows( L, 1 );
L1r = L1[ 1; 1 ];
L1i = L1[ 1; 2 ];

# Calculate the corresponding
# eigenvector and split it
# in its real and imaginary
# parts
V1 = matrixEigVect( A, L1 );
V1r = matrixGetColumns( V1, 1 );
V1i = matrixGetColumns( V1, 2 );

# We have A * V = V * diag( L )
# Thus the error for the first
# eigenvalue and eigenvector
# pair is
# A * V1 - V1 * L1

# Note that V1 and L1 have both
# real and imaginary parts so we
# need to do multiplication using
# (a+ib)(c+id)=(ac-bd)+(ad+cb)i

# Real portion of error
Er = A * V1r -
     ( V1r * L1r - V1i * L1i );

# Imaginary portion of error
Ei = A * V1i -
     ( V1i * L1r + V1r * L1i );

# Output
print( "A", A );
print( "L", L );
print( "L1", L1 );
print( "V1", V1 );

print( "Error" );

print ( "Er", Er );
print ( "Ei", Ei );