Here's a subroutine I wrote to recursively (!j) estimate Stock and Watson (2002) I(0) factors and Bai (2004) I(1) factors based on a Gauss code by Igor Masten where N>T.

Note: X is our matrix of time series created from g1_standstat1 - a group of standardized stationary variables.

Code: Select all

`'Stock and Watson Code:`

%start = "1980q1"

%endminusone = "2002q4-1"

!noiterations = 20

for !j = 1 to !noiterations

smpl {%start} {%endminusone}+!j

stom(g1_standstat1,X)

scalar T = @rows(X)

scalar N = @columns(X)

'maximum number of factors:

scalar k_max = 4

'Because N>T:

sym XX = X*@transpose(X)

scalar il = T-1-k_max+1

scalar iu = T-1

'Note: VA corresponds to VEctors, VE corresponds to VAlues

matrix va = @eigenvectors(XX)

vector ve = @eigenvalues(XX)

for !i = il to iu

vector factor!i = ((T-1)^(0.5))*@columnextract(va,!i)

matrix(T,k_max) factors

colplace(factors,factor!i,!i-il+1)

next

matrix loadings = @transpose(factors)*X/(T-1)

'For a normalization, multiply/divide by the element [kmax,1] of the loadings matrix.

scalar norm = loadings(k_max,1)

factors = factors*norm

loadings = loadings/norm

'Re-extract normalized factors from this matrix as series for future subroutines, in reverse order:

'i.e. the factor associated with the largest eigenvalue is F_1, second is F_2,...,F_k_max

mtos(factors,factorgroup1)

for !i = 1 to k_max

scalar temp = k_max-!i+1

rename ser0{temp} f_1_!i_!j

next

delete temp

delete factorgroup1

This is our subroutine for estimating I(1) factors, where g1_standnonstat1 is our group of nonstationary standardized variables and the subscript denotes the I(1) property.

Code: Select all

`'Bai (2004) Code`

stom(g1_standnonstat1,i_X)

scalar i_T = @rows(i_X)

scalar i_N = @columns(i_X)

scalar i_k_max = 4

sym i_XX = i_X*@transpose(i_X)

'Note the difference with the stationary routine here:

scalar i_il = T-k_max+1

scalar i_iu = T

matrix i_va = @eigenvectors(i_XX)

vector i_ve = @eigenvalues(i_XX)

for !i = i_il to i_iu

'Again, note the difference with the stationary routine here:

vector i_factor!i = (T)*@columnextract(i_va,!i)

matrix(i_T,i_k_max) i_factors

colplace(i_factors,i_factor!i,!i-i_il+1)

next

matrix i_loadings = @transpose(i_factors)*X/(i_T^2)

'normalize

scalar i_norm = i_loadings(i_k_max,1)

i_factors = i_factors*i_norm

i_loadings = i_loadings/i_norm

'Re-extract normalized factors from this matrix, in reverse order: i.e. the factor associated with the largest eigenvalue is if1, second is i_F_2,...,i_F_km_max

mtos(i_factors,i_factorsgroup1)

for !i = 1 to i_k_max

scalar temp = i_k_max-!i+1

rename ser0{temp} if1_!i_!j

next

delete temp

delete i_factorsgroup1

Best wishes