#
# Cut-and-Paste Code below into Window above and Run
#
require(dse)
require(matlab)
#
# WL20 BAU Model
#
# Measurement Matrix (Growth-LP), (P.Wheat.+LP-TEMP), (P.Wheat.+ P.OIL.- EF- Earths)
#
# N OIL QA GWP P.Wheat. P.Oil. TEMP CO2
#[1,] 0.2845 0.271 0.2821 0.2824 0.236 0.243 0.251 0.2837
#[2,] -0.0161 0.241 0.0985 -0.1413 0.411 0.066 -0.235 -0.0983
#[3,] -0.1112 -0.241 -0.0883 0.0213 0.461 0.712 0.190 -0.0201
# Carbon TotalFootprint Earths WorldGlobal LivingPlanet URBAN
#[1,] 0.2833 0.275 0.279 0.2819 -0.1839 0.2840
#[2,] 0.0925 0.182 0.148 -0.0926 0.7705 -0.0868
#[3,] -0.0955 -0.293 -0.231 -0.0943 -0.0576 -0.0707
#
#Fraction of Variance
#[1] 0.874 0.941 0.973 0.986 0.995 0.998 0.999 0.999 1.000 1.000 1.000
#[12] 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.01169571, -0.02519461, -0.09358484, 0.210941578,
-0.01360718, 0.96330129, -0.05250691, -0.003862528,
0.02303607, -0.03631204, 0.93554760 , 0.018115378,
0.000000000, 0.00000000, 0.0000000, 1.00000000
),byrow=TRUE,nrow=4,ncol=4)
#
# Creating a Steady State
#
# Eliminating (P.Wheat.+LP-TEMP) Feedbacks
#
# f[1,2] <- f[2,1] <- f[3,2] <- 0
#
#
h <- eye(3,4)
k <- f[1:4,1:3,drop=FALSE]
WL20 <- SS(F=f,H=h,K=k,z0=c( 0.210941578, -0.003862528, 0.018115378,
1.0000000),
output.names=c("W1","W2","W3"))
stability(SS(F=f[1:3,1:3,drop=FALSE],H=eye(3),Q=eye(3),R=eye(3)) )
shockDecomposition(toSSChol(WL20))
WL20.data <- simulate(WL20,sampleT=50,start=1950,noise=matrix(0,50,3),)
m <- l(WL20,WL20.data)
# tfplot(m)
AIC(m)
WL20.f <- forecast(m,horizon=50)
WL20.fx <- merge.forecast(WL20.f)
tfplot(WL20.f)
#
# MX_CA_US BAU Growth W Input
#
f <- matrix( c( 0.62800004, 0.01354979, 0.14299312, -0.1642985,
0.01480452, 0.71389792, -0.01375403, -0.1647911,
0.07563929, -0.04559007, 0.67184221, 0.1091528,
0.00000000, 0.00000000 , 0.00000000, 1.0000000
),nrow=4,ncol=4, byrow=TRUE)
g <- matrix( c( 0.08549373, 0.2251402, -0.1320120,
0.26038765 , 0.1120619, -0.1408007,
0.47578997 ,-0.2712762, 0.1180161,
0.00000000, 0.0000000, 0.0000000
),nrow=4,ncol=3, byrow=TRUE)
h <- eye(3,4)
k <- (f[,1:3,drop=FALSE])
MX_CA_US <- SS(F=f,H=h,G=g,K=k,z0=c( -0.6533331, -1.2629616, -1.5494593,1),
output.names=c("MX1","CA1", "US1"),
input.names=c("W1", "W2", "W3"))
print(MX_CA_US)
#is.SS(MX_CA_US)
stability(SS(F=f[1:3,1:3,drop=FALSE],H=eye(3),Q=eye(3),R=eye(3)) )
MX_CA_US.data <- simulate(MX_CA_US,sampleT=50,noise=matrix(0,50,3),
start=1950,input=WL20.fx)
seriesNames(outputData(MX_CA_US.data)) <- c("MX1","CA1", "US1")
tfplot(MX_CA_US.data)
shockDecomposition(toSSChol(MX_CA_US),shock=-1)