demo/dyn-mcmcpack.R

In dyn: Time Series Regression

# based on code by Delali Accolley developed for:
# Time Series Analysis: Some Applications, Dec 30, 2008.
# http://www.scribd.com/people/view/5307

# Money Supply
m <- c(2.14519640611418, 2.14983469671578, 2.14921911265538, 
2.14612803567824, 2.14550717140966, 2.14488541828714, 2.14983469671578, 
2.14829409743475, 2.15198239545747, 2.15503222879097, 2.15775888604686, 
2.16196661636407, 2.16435285578444, 2.16613397030511, 2.16524432612531, 
2.16967443405881, 2.17376882313665, 2.17724783625562, 2.18184358794477, 
2.18554215485438, 2.18892848376085, 2.19200959265367, 2.20057692675485, 
2.20493352235414, 2.20817252666712, 2.21005084987514, 2.21722065564452, 
2.22479195649268, 2.23172438332852, 2.23451728351269, 2.23552844690755, 
2.23552844690755, 2.24254142829838, 2.24797326636181, 2.25695815256093, 
2.26316246496222, 2.26834391395106, 2.27737997466725, 2.28488171465545, 
2.29534714833362, 2.30102999566398, 2.30384377488865, 2.30556631351530, 
2.30941722577814, 2.31323429169472, 2.31722734917642, 2.32592595577147, 
2.33122478102073, 2.34004731766139, 2.34927752746796, 2.35506820634885, 
2.35850591149024, 2.36977228859696, 2.37401474029191, 2.38596357060070, 
2.39654803798713, 2.40088321554836, 2.40942586867144, 2.41145134213794, 
2.41979058610636, 2.42602301568988, 2.42894429003557, 2.43296929087441, 
2.43806745045349, 2.44153803870216, 2.45086469237977, 2.45591024038274, 
2.45803319249651, 2.46642272243379, 2.47158505418519, 2.47654180902743, 
2.48600518636224, 2.49679131570004, 2.50351831272407, 2.51121470113639, 
2.51969676715985, 2.52762990087134, 2.54020429984206, 2.54678935163126, 
2.55303301620244, 2.55930801090701, 2.57217431361306, 2.57898284270279, 
2.58183592405765, 2.58927922123597, 2.59006123080374, 2.60724050383174, 
2.61119206086843, 2.62221402296630, 2.62859325585126, 2.63042787502502, 
2.64018319192134, 2.64610952197885, 2.65079303965193, 2.65982115805570, 
2.67651071028255, 2.69072754387037, 2.70243053644553, 2.70994801651076, 
2.71717102683231, 2.72542155007426, 2.73255457985143, 2.73487980279263, 
2.74162425750381, 2.75327657018442, 2.76514679010803, 2.78053332531640, 
2.79225157190326, 2.80174661921946, 2.81993856935540, 2.83720952780121, 
2.86015826131828, 2.86557770741993, 2.87110570098559, 2.87366929270679, 
2.87517705981470, 2.88184096832493, 2.89114703044877, 2.89414984676792, 
2.89580915016913, 2.89376176205794, 2.88846031803539, 2.8926510338773, 
2.89921841785137, 2.90395770852317, 2.90789483541628, 2.91386681189624, 
2.91629599456313, 2.92360664301746, 2.9328794578238, 2.9396190789567, 
2.95279244304409, 2.97160052023006, 2.97968492382703, 2.99480089929460, 
3.01068149313144, 3.01640650087112, 3.03140846425162, 3.04304774280412, 
3.05292468370773, 3.05701912432277, 3.05884341466876, 3.06141477827707, 
3.06096211309383, 3.05941193743866, 3.05846398560225, 3.05762807295557, 
3.05207803048417, 3.05018634967536, 3.04735276075393, 3.03985017774966, 
3.03394620299036, 3.03031630598859, 3.02767571590489, 3.02840856511547, 
3.03039730085676, 3.03237697120994, 3.0321350468799, 3.03338354117312, 
3.03961238189672, 3.04040452891416, 3.04119523369681, 3.04000863601354, 
3.05018634967536, 3.04469651711019, 3.04277233749767, 3.04103720786703, 
3.03638932866569, 3.04536210265335, 3.05176974689933, 3.08066255639418, 
3.07258073264895, 3.07678603350808, 3.07649480500972, 3.07729521393209, 
3.08625302381716, 3.09335163201555, 3.10748131891124, 3.11270552103089, 
3.11597642945482, 3.12404751911003, 3.12804367382646, 3.13382621407639, 
3.13871310987601, 3.13776542005734, 3.13987908640124, 3.13896547828983, 
3.13814474417949, 3.1414811292708, 3.13842902000600, 3.13385812520333, 
3.13560963602868, 3.1365303235313, 3.13529170447575, 3.13538710838239, 
3.13548249133571, 3.13820793266816, 3.14160653011825)

# Gross Domestic Product
gdp <- c(3.37890384090191, 3.39015967351696, 3.38989279010299, 
3.39138684913375, 3.40092784656566, 3.39876008920108, 3.3994308587464, 
3.39378167908042, 3.39638381539554, 3.40449968813707, 3.41145400372149, 
3.4202204669763, 3.42797190383071, 3.43270894117274, 3.43669995864986, 
3.43775627244322, 3.44338720522461, 3.44883391975793, 3.45688964600751, 
3.46030022342821, 3.46990168182446, 3.47488435713166, 3.480827066323, 
3.48196552361768, 3.49250952954897, 3.49830601521825, 3.50706672648902, 
3.51745526225994, 3.52795477974597, 3.52941509669321, 3.53229660223629, 
3.53575970117355, 3.53957295156362, 3.53963324001493, 3.54302597193560, 
3.54630382645788, 3.55516088852466, 3.56246826031581, 3.56542206550865, 
3.56723494995109, 3.57404765011263, 3.57530536900184, 3.57799015802287, 
3.57589525815796, 3.57517903361618, 3.57598392316098, 3.57983700320023, 
3.57516085343014, 3.58702814961348, 3.58949902807404, 3.59290482679588, 
3.59415022081285, 3.60182514332519, 3.61193208097921, 3.61607608380287, 
3.62311406263558, 3.63399789876006, 3.63900638749519, 3.63669516595899, 
3.64082136150782, 3.63703554559244, 3.63827470655617, 3.6340590351007, 
3.63236003553766, 3.62711730105762, 3.63028553004575, 3.63756993685921, 
3.64324745713485, 3.65291175275886, 3.65613381518738, 3.65819414921528, 
3.661310463251, 3.66651131603438, 3.67495278091234, 3.68267314451279, 
3.68262719974132, 3.68402867871117, 3.70080027423614, 3.70507295407361, 
3.71074471499261, 3.71159062045203, 3.71200116361905, 3.71512078978917, 
3.71639450720212, 3.71778208606288, 3.70892824825998, 3.70819629117332, 
3.71618401341002, 3.72488790086750, 3.7214953493211, 3.72670680981111, 
3.72126570730132, 3.71409006177185, 3.71640269185504, 3.7147742873117, 
3.71515696446622, 3.72048220117837, 3.73015786224485, 3.73865575592748, 
3.74744923480017, 3.75586437919688, 3.76327048676597, 3.76747542182771, 
3.7710225866597, 3.77501903753436, 3.77871412737249, 3.78545156010774, 
3.78877606475967, 3.79290455793201, 3.79463126861595, 3.79876829073298, 
3.80095056584114, 3.8037986386069, 3.80855196495661, 3.81247348262485, 
3.81999496917294, 3.82210639725106, 3.82759579522881, 3.82990240971491, 
3.83560026578289, 3.83998991188877, 3.84282601971825, 3.84591237405979, 
3.84700987753756, 3.85200054395039, 3.85310643693097, 3.85313225606482, 
3.84984335552428, 3.84762647722758, 3.85042743626793, 3.85252181617425, 
3.85455900268734, 3.8590316961685, 3.86320016532147, 3.86743813475053, 
3.87219860718782, 3.87272126540949, 3.87491957807085, 3.87713908978218, 
3.88294580194789, 3.88733904638370, 3.89296997868462, 3.89539088656351, 
3.90045674168868, 3.90166014684618, 3.90243682014472, 3.90596207489974, 
3.90912601325841, 3.91217662491938, 3.91923838972054, 3.92286403948605, 
3.92791019872975, 3.93125527181612, 3.93780931723206, 3.94318587698685, 
3.94637230173611, 3.95115171225313, 3.95401453446006, 3.95898614490071, 
3.96553349517956, 3.96920799549647, 3.97278494082567, 3.97782702104881, 
3.98547606390483, 3.98657563478803, 3.99334286067809, 3.99284603331591, 
3.99509737189303, 3.99456083155732, 3.9958904721239, 3.99436332657822, 
3.99607670797315, 3.99901344848855, 4.00136835707199, 4.00392024283679, 
4.00413660082482, 4.00543742651056, 4.00914082023063, 4.01697677473657, 
4.01981923198026, 4.02298742497743, 4.0267083073446, 4.03054851765607, 
4.03327955375252, 4.03646343237022, 4.03925921694932, 4.04336179309352, 
4.04478029595905, 4.04988554653412, 4.05275712606633, 4.05361991176023, 
4.05524013376498, 4.05529792736715, 4.06037265619674, 4.06542249799062, 
4.06523285288133, 4.06617612040347, 4.06913320016068)

library(dyn) # also pulls in zoo
library(MCMCpack) # also pulls in coda

m. <- zooreg(m, start = 1959, freq = 4)
gdp. <- zooreg(gdp, start = 1959, freq = 4)

L <- function(x, k) lag(x, -k)

summary(post <- dyn$MCMCregress(m. ~ L(m., 1:3) + L(gdp., 1:3)))
varnames(post) <- c("Intercept", "M[t-1]", "M[t-2]", "M[t-3]", 
	"GDP[t-1]", "GDP[t-2]", "GDP[t-3]", "sigma^2")

# posterior and diagnostics
summary(post)
raftery.diag(post)
geweke.diag(post)

# Density plot
par(mfrow = c(4, 2))
for(v in varnames(post)) 
	densplot(post[, v], col = "red", main = parse(text = v))
par(mfrow = c(1, 1))

# Trace plot 
par(mfrow = c(4, 2))
for(v in varnames(post)) 
	traceplot(post[, v], col = "blue", main = parse(text = v))
par(mfrow = c(1, 1))