#
# World1 Forrester (Cut and paste code into window above and Run (Cmd-Enter)
#
AIC <- function(model) {informationTestsCalculations(model)[3]}
require(dse)
require(matlab)
#
# Measurement Matrix
# NR QL N K POL
#[1,] -0.46659 -0.399 0.440 0.464 0.4619
#[2,] -0.00684 0.841 0.494 0.220 0.0279
#[3,] -0.36164 0.354 -0.596 -0.105 0.6143
#
# Fraction of Variance
# 0.909 0.984 1.000 1.000 1.000
#
f <- matrix( c(0, 0, 0, -0.93989620, 0.09577918, -0.03382348, 0.000000000,
0, 0, 0, -0.06989283, -1.13162828, 0.07219068, 0.000000000,
0, 0, 0, -0.08617730, -0.19537459, -0.90189005, 0.000000000,
1, 0, 0, 1.93850140, -0.09551917, 0.03364973, 0.002188603,
0, 1, 0, 0.07149041, 2.13136386, -0.07368290, -0.003111204,
0, 0, 1, 0.08816385, 0.19670294, 1.89730057, -0.003987851,
0, 0, 0, 0.00000000, 0.00000000, 0.00000000, 1.000000000
),byrow=TRUE,nrow=7,ncol=7)
k <- (f[1:7,3:5,drop=FALSE])
h <-matrix( c(0, 0, 0, 1, 0, 0, 0,
0 , 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 1 , 0
),byrow=TRUE,nrow=3,ncol=7)
W1 <- SS(F=f,H=h,K=k,z0=c(0, 0, 0, 0.002188603, -0.003111204,-0.003987851,1),
output.names=c("W1","W2","W3"))
print(W1)
is.SS(W1)
stability(W1)
# tfplot(simulate(W1,sampleT=100))
W1.data <- simulate(W1,sampleT=100,noise=matrix(0,100,3),start=1900)
W1.f <- forecast(l(W1,W1.data),horizon=100)
tfplot(W1.f)
shockDecomposition(toSSChol(W1))
AIC(l(W1,W1.data))