[R] ts model challenge (transfer function)

[Felix Andrews]

Dear useRs,

I am trying to model a time series with a transfer function. I think
it can be put into the ARMA framework, and estimated with the 'arima'
function (and others have made similar comments on this list). I have
tried to do that, but the results have so far been disappointing.
Maybe I am trying to make 'arima' do something it can't...

The data are time series of river flow 'x' (output) and rainfall 'u'
(exogenous input).

x <- c(0.645, 0.61, 0.594, 0.579, 0.558, 0.62, 3.977, 2.824, 0.985,
0.764, 0.691, 0.659, 0.633, 0.619, 0.613, 0.708, 0.692, 0.617,
0.599, 0.586, 0.574, 0.566, 0.565, 0.558, 0.588, 0.6, 0.566,
0.887, 0.77, 0.692, 0.846, 0.829, 0.65, 0.613, 0.59, 1.256, 0.85,
0.673, 0.644, 0.649, 0.694, 0.663, 0.629, 0.736, 0.692, 0.611,
0.595, 0.604, 0.626, 0.618, 1.501, 1.255, 0.856, 1.742, 1.575,
1.003, 0.912, 0.972, 1.229, 1.812, 4.067, 6.921, 5.866, 3.26,
2.129, 1.744, 1.566, 2.445, 2.362, 1.77, 1.594, 1.713, 1.784,
1.56, 1.41, 1.839, 3.008, 4.697, 10.853, 7.325, 3.853, 2.96,
2.54, 2.315, 3.488, 4.414, 2.782, 3.498, 6.411, 7.902, 8.18,
7.468, 4.437, 3.663, 3.255, 2.972, 2.757, 2.575, 2.443, 2.316,
2.226, 2.223, 2.166, 2.154, 3.346, 4.047, 4.059, 3.134, 2.628,
2.439, 2.315, 2.238, 2.152, 2.116, 2.31, 3.864, 5.727, 5.356,
4.582, 5.626, 5.383, 4.368, 4.063, 3.963, 3.604, 3.591, 3.25,
3.122, 3.585, 3.488, 3.063, 2.865, 2.672, 2.946, 3.415, 3.237,
2.715, 2.831, 2.59, 2.441, 2.368, 2.285, 2.272, 2.633, 2.332,
2.199, 2.135, 2.063, 2.019, 1.998, 1.967, 1.945, 1.986, 2.004,
2.047, 2.735, 2.101, 2.038, 1.978, 1.94, 1.983, 2.442, 2.239,
2.081, 3.18, 3.039, 2.134, 1.969, 1.857, 1.794, 1.771, 1.703,
1.687, 1.658, 1.64, 1.64, 1.589, 1.546, 1.52, 1.515, 1.505, 1.505,
1.484, 1.477, 1.488, 1.709, 1.505, 1.522, 2.309, 1.836, 1.528,
1.453, 1.417, 1.373, 1.331, 1.307, 1.308, 1.338, 1.29, 1.292,
1.291, 1.466, 1.438, 1.358, 1.31, 1.25, 1.224, 1.192, 1.159,
1.125, 1.107, 1.093, 1.113, 1.114, 1.132, 1.125, 1.077, 1.139,
1.156, 1.082, 1.045, 1.021, 1.121, 1.735, 1.326, 1.126, 1.049,
1.111, 1.08, 1.027, 1, 0.978, 0.947, 0.946, 0.97, 0.936, 0.92,
0.904, 0.878, 0.887, 1.024, 0.994, 1.007, 0.915, 0.894, 0.834,
0.82, 0.799, 0.812, 0.814, 0.786)

u <- c(0, 0, 0, 0, 0.003, 0, 32.551, 3.369, 0, 0.634, 0, 0.129, 0,
0, 0.052, 0.959, 0.634, 0, 0, 0, 0, 0, 0, 0.011, 0.266, 0.116,
0.077, 1.976, 0.288, 0.63, 0.514, 0, 0, 0, 4.062, 0.012, 0, 0.143,
0.291, 0.53, 0.018, 0.14, 0.99, 0.143, 0.199, 0, 0.202, 0.167,
0.069, 3.413, 0.731, 1.892, 2.793, 3.787, 0.658, 0.717, 1.63,
2.375, 4.119, 8.937, 17.177, 11.645, 17.018, 0.54, 0.021, 0.354,
0, 9.792, 2.851, 1.541, 2.498, 2.042, 3.078, 0, 0.119, 5.864,
10.555, 17.609, 28.071, 21.461, 0, 0, 0.027, 0.776, 12.466, 6.447,
8.693, 9.736, 27.098, 26.849, 5.95, 0.058, 0, 0, 0.068, 0.015,
0.07, 0.308, 0.066, 0.013, 0.721, 0.606, 0.792, 9.646, 26.848,
2.4, 0.112, 1.749, 0, 0.183, 0.065, 0, 0.256, 0.932, 11.855,
18.217, 19.209, 7.409, 17.594, 16.336, 2.439, 3.926, 3.677, 1.26,
2.85, 0.833, 0.094, 4.78, 5.436, 0, 0, 0.183, 5.006, 3.348, 0.362,
2.042, 1.653, 0.751, 0, 0.023, 0, 0.43, 4.654, 0.012, 0, 0, 0,
0.117, 0, 0.087, 0.005, 0.701, 0, 0.506, 1.257, 2.697, 1.028,
0, 0.078, 0, 0.39, 4.624, 0.186, 1.925, 6.711, 5.125, 0.112,
0.121, 0.006, 0, 0.344, 0, 0.003, 0, 0.572, 0.007, 0, 0, 0, 0,
0.007, 0.03, 0, 0.006, 1.439, 0, 1.007, 3.312, 0.026, 0.094,
0, 0.093, 0, 0, 0, 0.45, 0.011, 0, 0, 0, 0.777, 0.077, 0.184,
0.1, 0, 0.038, 0, 0, 0, 0, 0.009, 0, 0.009, 0.002, 0.087, 0,
0, 0.21, 0.056, 0.046, 0, 0, 0.578, 1.17, 0.426, 0, 0.318, 0.124,
0.096, 0, 0, 0, 0.012, 0, 0, 0.002, 0, 0, 0, 0.053, 0.178, 0.068,
0.126, 0, 0.001, 0.001, 0, 0.001, 0, 0, 0)

matplot(cbind(flow,rain), type="lh")
# Note that 'u' is highly skewed.

# The system can be modelled by this transfer function:
# x[i] <- a_1*x[i-1] + a_2*x[i-2] + b_0*u[i] + b_1*u[i-1]

# I happen to know that a good model is
a_1 <- 1.6545
a_2 <- -0.6580
b_0 <- 0.1149
b_1 <- -0.1115
# (This was estimated by a proprietary program using
#  "Simple Refined Instrumental Variable" algorithm).

# The transfer function differs from the model used by 'arima', as
# there is a scale factor applied to the instantaneous input u[i].
# So need to scale the input by b_0 and also scale b_1:
arOK <- c(a_1, a_2)
maOK <- b_1 / b_0
scaleOK <- b_0

sim.good <- arima.sim(model=list(ar=arOK, ma=maOK), innov=u*scaleOK,
n=length(x))
lines(sim.good, col="blue")
# good fit to data

# Now I try to estimate the same model using 'arima'...
# Can I use xreg=u to estimate the scale factor b_0 ?
cal <- arima(x=x, order=c(2,0,1), xreg=u, include.mean=F)
round(coef(cal),digits=4)
#    ar1    ar2    ma1      u
# 0.4116 0.4965 0.6992 0.0644

model.new <- list(ar=coef(cal)[1:2], ma=coef(cal)[3])
scale.new <- coef(cal)[["u"]]
sim.new <- arima.sim(model=model.new , innov=u*scale.new , n=length(x))
lines(sim.new , col="green")
# poor fit to data

This result is disappointing, considering how good the other parameter
set is. I think there must be something I am not understanding about
'arima' -- probably 'xreg'.

Many thanks for any help.

Felix

-- 
Felix Andrews / 安福立
PhD candidate, The Fenner School of Environment and Society
The Australian National University (Building 48A), ACT 0200
Beijing Bag, Locked Bag 40, Kingston ACT 2604
http://www.neurofractal.org/felix/
voice:+86_1051404394 (in China)
mobile:+86_13522529265 (in China)
mobile:+61_410400963 (in Australia)
xmpp:foolish.android at gmail.com
3358 543D AAC6 22C2 D336  80D9 360B 72DD 3E4C F5D8