Snippets
Model
Smoothing
Simple and exponential
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# This constructs a simple
# weighted and exponentially
# weighted moving average
# series for US CPI.
# Two exponentially weighted
# moving averages are
# constructed, one using
# a weight parameter and
# one with a length
# parameter.
CPI = dataFred( "CPIAUCSL" );
# Percentage change (mom)
pCPI = p( CPI );
# Smooth
s1 = sma( 12, pCPI );
s2 = ewmaLength( 12, pCPI );
s3 = ewma( 0.99, pCPI );
# Throw them all together in
# an observation matrix
om = matrixObs( CPI, s1, s2, s3 );
# Print last few rows
c = cols( om );
r = rows( om );
print( om[ r-5:r ; 1:c ] );
# Use last smoothed value
# to generate forecastss
l = lastDate( pCPI );
print( "Final value",
observe( l, s1[ l ] ),
observe( l, s2[ l ] ),
observe( l, s3[ l ] ) );
p1 = CPI[ l ] * ( 1 + s1[ l ] );
p2 = CPI[ l ] * ( 1 + s2[ l ] );
p3 = CPI[ l ] * ( 1 + s3[ l ] );
print( "Forecast",
p1, p2, p3 );
print( "YOY%",
100 * ( p1 / CPI[ l - 11 ] - 1 ),
100 * ( p2 / CPI[ l - 11 ] - 1 ),
100 * ( p3 / CPI[ l - 11 ] - 1 ) );
OLS
Fit an OLS model
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# This constructs a simple oil
# price model - best run in
# monthly frequency. It does
# some calculations with the
# model and also uses it to
# forecast out of sample.
# Basic format is
# ols( "Name", y, 1, x1, x2 ... )
# The model is
# dLn( p ) = a0 +
# a1 * dLn( p(-1) ) +
# a2 * dLn( x(-1) ) +
# a3 * dLn( d(-1) )
# Where
# p is the price of oil
# x is the denominating currency
# d is demand, here proxied by
# industrial production
# Lower accuracy to speed up
# calculations a bit
# (default accuracy is 50 digits)
accuracy( 15 );
# Data
oilPrice = dataFred( "WTISPLC" );
usdIndex = dataFred( "TWEXBMTH" );
indProd = dataFred( "INDPRO" );
# Precalculate as we'll use
# these a lot
Y = dLn( oilPrice );
X1 = dLn( usdIndex );
X2 = dLn( indProd );
# Fit the model from 1995 onwards
# Use arbitrary large ending index
i0 = dateIndex( 1995, 1, 1 );
i1 = dateIndex( 2100, 1, 1 );
# Store the model as "oilModel"
ols( "oilModel",
# This limits the sample
limit( i0, i1, Y ),
1, # Constant
Y(-1),
X1(-1),
X2(-1) );
# Now we use the saved results to
# calculate the variance inflation
# factor for each regressor
n = dataSetValue( "oilModel",
"nObs" );
# First VIF
# Hide output with '-' prefix
ols( "-a1",
limit( i0, i0 + n - 1, Y(-1) ),
1,
X1(-1),
X2(-1) );
R21 = dataSetValue( "a1",
"nRSq" );
# Second VIF
# Hide output with '-' prefix
ols( "-a2",
limit( i0, i0 + n - 1, X1(-1) ),
1,
Y(-1),
X2(-1) );
R22 = dataSetValue( "a2",
"nRSq" );
# Third VIF
# Hide output with '-' prefix
ols( "-a3",
limit( i0, i0 + n - 1, X2(-1) ),
1,
Y(-1),
X1(-1) );
R23 = dataSetValue( "a3",
"nRSq" );
print( "Variance inflation factors" );
print( "VIF1", 1 / ( 1 - R21 ) );
print( "VIF2", 1 / ( 1 - R22 ) );
print( "VIF3", 1 / ( 1 - R23 ) );
# Now use the model to project
l = lastDate( oilPrice );
yhat = dataSetValue(
"oilModel",
"projected" );
print( limit( l - 11,
l + 1,
yhat ) );
# The projected value is in
# change in logs, transform
# back
f = oilPrice[ l ] *
exp( yhat[ l + 1 ] );
print( "History",
limit( l - 12, l,
oilPrice ) );
print( "Forecast",
observe( l + 1, f ) );
RWA
Risk weighted assets
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Calculate capital
# requirement and risk
# weighted assets for a
# hypothetical bank
# See “An explanatory note
# on the Basel II IRB risk
# weight functions” July
# 2005 (BIS)
# Three exposures
# 1. Residential mortgages
# 2. Revolving retail
# 3. Other retail
# Exposure at default
# (Arbitrary)
EAD = matrix( 3, 1,
50000000, # 50m
25000000, # 25m
15000000 ); # 15m
# Probability of default
# (Arbitrary)
PD = matrix( 3, 1,
0.15,
0.25,
0.30 );
# Loss given default
# (Arbitrary)
LGD = matrix( 3, 1,
0.35,
0.65,
0.85 );
# Maturity, in years
# (Arbitrary)
M = matrix( 3, 1,
5.5,
0.8,
0.3 );
# Smoothed maturity
# adjustment
bPD =
matrixElPow(
matrixC( 0.11852, 3 ) -
0.05478 * ln( PD ),
matrixC( 2, 3 ) );
# Asset correlations
R = matrix( 3, 1,
0.15,
0.04,
0.03 *
( 1 - exp
( -35 * PD[ 3;1 ] ) ) /
( 1 - exp( -35 ) ) +
0.16 * ( 1 - ( 1 - exp
( -35 * PD[ 3;1 ] ) ) /
( 1 - exp( -35 ) ) ) );
# Capital requirement
K = matrix( 3, 1 );
i = 1;
while( i <= 3,
K[ i; 1 ] = ( LGD[ i;1 ] *
normPdf
( ( 1 - R[ i;1 ] ) ^ -0.5 *
normQnt( PD[ i;1 ] ) +
( R[ i;1 ] /
( 1 - R[ i;1 ] ) ) ^ 0.5
* normQnt( 0.999 ) ) -
PD[ i;1 ] * LGD[ i;1 ] ) /
( 1 - 1.5 * bPD[ i;1 ] ) *
( 1 + ( M[ i;1 ] - 2.5 ) *
bPD[ i;1 ] ),
i = i + 1 );
# Risk weighted asset
RWA = 12.5 * matrixElMul( K, EAD );
print(
"Capital requirement %",
K );
NomK = matrixElMul( K, EAD );
print(
"Nominal capital requirement",
NomK, matrixSum( NomK ) );
print( "RWA", RWA,
matrixSum( RWA ) );
ARIMA
Fit an ARIMA model
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# This constructs an ARIMA
# model to monthly US CPI. It
# also produces some forecasts
# using the model and then
# prepare a report with those.
# Basic format is
#
# arima( "Name", y, 1, x1, x2 ...
#
# to specify the LHS and any RHS
# variables followed by a matrix
# with up to 7 entries
# ( p, d, q, P, D, Q, S ).
# Since these are given on a single
# row, for e.g. ARIMA( 1, 1, 1 ) use
# matrix( 1, 3, 1, 1, 1 )
# as this parameter. This is followed
# by the forecast horizon and then
# something like "-exp" to transform
# the forecast back to levels if it was
# estimated in logs for example.
# The model is
# Data
CPI = dataFred( "CPIAUCSL" );
y = ln( CPI );
# Print the acf and pacf
bars( acf( d( y ) ) );
bars( pacf( d( y ) ) );
# Fit the model
# Note that here we use a
# ARIMA( 2, 1, 2 ) model.
# Note that we will store it to
# "cpi" so remove the report
# beforehand just in case it
# already exists. Also note that
# we backtransform the forecasts
# to CPI levels.
reportDelete( "cpi" );
arima( "#cpi", y, 1, matrix( 1, 3, 2, 1, 2 ), "-exp" );
ARIMA
US GDP ARIMA
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# This fits an ARIMA
# model to quarterly
# US CPI. Then it
# creates a report
# containing the
# model and some
# forecasts made by it.
# This assumes that the
# series GDPC1 is available
# as a date series
y = ln( GDPC1 );
# Print the acf and pacf
bars( acf( d( y ) ) );
bars( pacf( d( y ) ) );
# Fit the model
arima( "#GDP", y, 1, matrix( 1, 3, 1, 1, 1 ), "-exp" );