require(dse)

require(matlab)

#

# USL20W BAU Model

#

# Measurement Matrix 

#           Q      N  HOURS XREAL      X      L

#[1,]  0.4124  0.403  0.413 0.410  0.399  0.412

#[2,] -0.1464 -0.580 -0.122 0.310  0.711 -0.160

#[3,] -0.0825 -0.546  0.137 0.358 -0.466  0.575

#

# Fraction of Variance 

#[1] 0.977 0.999 1.000 1.000 1.000 1.000

#

merge.forecast <- function (fx,n=1) {

#

# Merges a forecast with the outputdata

x <- splice(fx$pred,fx$forecast[[n]])

colnames(x) <- seriesNames(fx$data$output)

return(x)

}

AIC <- function(model) {informationTestsCalculations(model)[3]}

f <- matrix( c(   1.026261160, -0.04873236, 0.04995449,  0.166789178,

                        0.023195483,  1.01292512 ,0.06694359,  0.004293523,

                       -0.002685722, -0.04074830, 1.04458649, -0.006797188,

                       0.000000000, 0.00000000,  0.0000000,  1.00000000

),byrow=TRUE,nrow=4,ncol=4)

h <- eye(3,4)

k <- f[1:4,1:3,drop=FALSE]

USL20 <- SS(F=f,H=h,K=k,z0=c( 0.166789178,  0.004293523, -0.006797188,  

1.0000000),

              output.names=c("US1","US2","US3"))

stability(USL20)

shockDecomposition(toSSChol(USL20))

USL20.data <- simulate(USL20,sampleT=50,start=1950)

m <- l(USL20,USL20.data)

# tfplot(m)

AIC(m)

USL20.f <- forecast(m,horizon=50)

USL20.fx <- merge.forecast(USL20.f)

tfplot(USL20.f)

#

# USL20c Growth Controlled

# by stabilizing growth rates

#

# f[1,1] <- 1.02

# f[2,2] <- 0.97

# f[3,3] <- 0.97

#

USL20c <- SS(F=f,H=h,K=k,z0=c( 0.166789178,  0.004293523, -0.006797188,  

1.0000000),

              output.names=c("US1","US2","US3"))

stability(USL20c)

shockDecomposition(toSSChol(USL20c))

USL20c.data <- simulate(USL20c,sampleT=50,start=1950)

m <- l(USL20c,USL20c.data)

# tfplot(m)

AIC(m)

USL20c.f <- forecast(m,horizon=50)

USL20c.fx <- merge.forecast(USL20c.f)

tfplot(USL20c.f)