Skip to content

Instantly share code, notes, and snippets.

@alexchinco
Created February 19, 2015 13:39
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save alexchinco/a34f3d0f34a1820e5ba5 to your computer and use it in GitHub Desktop.
Save alexchinco/a34f3d0f34a1820e5ba5 to your computer and use it in GitHub Desktop.
Data and code for post on impulse-response functions for VARs.
DATE VALUE
1947-01-01 .
1947-04-01 1.42857
1947-07-01 2.18083
1947-10-01 2.83237
1948-01-01 2.11874
1948-04-01 1.59207
1948-07-01 1.68382
1948-10-01 -0.91815
1949-01-01 -0.95147
1949-04-01 -0.10859
1949-07-01 -0.83623
1949-10-01 -0.24033
1950-01-01 -0.30854
1950-04-01 0.76313
1950-07-01 1.83448
1950-10-01 2.02454
1951-01-01 4.06593
1951-04-01 0.97288
1951-07-01 -0.05396
1951-10-01 1.48074
1952-01-01 0.37998
1952-04-01 0.26498
1952-07-01 0.67958
1952-10-01 0.11250
1953-01-01 -0.28842
1953-04-01 0.37566
1953-07-01 0.46033
1953-10-01 0.17509
1954-01-01 0.23429
1954-04-01 -0.15954
1954-07-01 -0.26013
1954-10-01 -0.30924
1955-01-01 0.13454
1955-04-01 -0.13436
1955-07-01 0.07475
1955-10-01 0.29876
1956-01-01 0.01117
1956-04-01 0.65897
1956-07-01 1.03562
1956-10-01 0.85295
1957-01-01 0.82396
1957-04-01 0.84962
1957-07-01 0.89244
1957-10-01 0.48473
1958-01-01 1.18662
1958-04-01 0.67161
1958-07-01 -0.05876
1958-10-01 0.10376
1959-01-01 0.17275
1959-04-01 0.17246
1959-07-01 0.51648
1959-10-01 0.60631
1960-01-01 0.09193
1960-04-01 0.59870
1960-07-01 0.05748
1960-10-01 0.64211
1961-01-01 0.20148
1961-04-01 -0.03351
1961-07-01 0.39222
1961-10-01 0.14359
1962-01-01 0.39013
1962-04-01 0.37533
1962-07-01 0.28789
1962-10-01 0.24087
1963-01-01 0.31929
1963-04-01 0.18375
1963-07-01 0.61245
1963-10-01 0.27018
1964-01-01 0.41230
1964-04-01 0.16166
1964-07-01 0.22595
1964-10-01 0.46055
1965-01-01 0.31097
1965-04-01 0.63918
1965-07-01 0.29533
1965-10-01 0.52877
1966-01-01 0.93543
1966-04-01 0.90492
1966-07-01 0.86588
1966-10-01 0.81553
1967-01-01 0.25545
1967-04-01 0.60667
1967-07-01 1.00401
1967-10-01 1.09552
1968-01-01 0.98326
1968-04-01 0.97368
1968-07-01 1.35233
1968-10-01 1.23714
1969-01-01 1.22485
1969-04-01 1.57805
1969-07-01 1.37238
1969-10-01 1.53521
1970-01-01 1.60000
1970-04-01 1.39895
1970-07-01 1.03538
1970-10-01 1.45262
1971-01-01 0.84091
1971-04-01 0.91904
1971-07-01 0.99256
1971-10-01 0.73710
1972-01-01 0.81220
1972-04-01 0.64597
1972-07-01 0.80048
1972-10-01 1.03498
1973-01-01 1.57198
1973-04-01 2.09142
1973-07-01 1.97346
1973-10-01 2.52902
1974-01-01 2.97607
1974-04-01 2.67865
1974-07-01 2.81261
1974-10-01 3.07212
1975-01-01 2.13729
1975-04-01 1.20418
1975-07-01 2.00564
1975-10-01 1.84274
1976-01-01 1.14535
1976-04-01 0.89445
1976-07-01 1.59574
1976-10-01 1.45375
1977-01-01 1.83545
1977-04-01 1.74493
1977-07-01 1.38462
1977-10-01 1.47379
1978-01-01 1.72027
1978-04-01 2.27500
1978-07-01 2.32677
1978-10-01 2.32389
1979-01-01 2.51852
1979-04-01 3.17919
1979-07-01 3.22129
1979-10-01 3.16554
1980-01-01 3.94566
1980-04-01 3.37454
1980-07-01 1.87638
1980-10-01 2.80418
1981-01-01 2.76508
1981-04-01 2.08568
1981-07-01 2.78499
1981-10-01 1.62572
1982-01-01 0.88837
1982-04-01 1.44503
1982-07-01 1.73601
1982-10-01 0.30727
1983-01-01 0.06841
1983-04-01 1.15612
1983-07-01 0.97546
1983-10-01 0.99900
1984-01-01 1.41741
1984-04-01 0.94311
1984-07-01 0.86957
1984-10-01 0.86207
1985-01-01 0.91833
1985-04-01 0.90903
1985-07-01 0.62201
1985-10-01 1.01946
1986-01-01 0.52018
1986-04-01 -0.48737
1986-07-01 0.61174
1986-10-01 0.69918
1987-01-01 1.20670
1987-04-01 1.13327
1987-07-01 1.06132
1987-10-01 0.93290
1988-01-01 0.78035
1988-04-01 1.14769
1988-07-01 1.21888
1988-10-01 1.09244
1989-01-01 1.13633
1989-04-01 1.61589
1989-07-01 0.78215
1989-10-01 1.01685
1990-01-01 1.72086
1990-04-01 0.98959
1990-07-01 1.72699
1990-10-01 1.69843
1991-01-01 0.74757
1991-04-01 0.59362
1991-07-01 0.76198
1991-10-01 0.82943
1992-01-01 0.67812
1992-04-01 0.76875
1992-07-01 0.76360
1992-10-01 0.87571
1993-01-01 0.72800
1993-04-01 0.72204
1993-07-01 0.46287
1993-10-01 0.82892
1994-01-01 0.50217
1994-04-01 0.56783
1994-07-01 0.92657
1994-10-01 0.58227
1995-01-01 0.73447
1995-04-01 0.81728
1995-07-01 0.50427
1995-10-01 0.54492
1996-01-01 0.88939
1996-04-01 0.85963
1996-07-01 0.57545
1996-10-01 0.86904
1997-01-01 0.60882
1997-04-01 0.22990
1997-07-01 0.50000
1997-10-01 0.53918
1998-01-01 0.20598
1998-04-01 0.32901
1998-07-01 0.51313
1998-10-01 0.46888
1999-01-01 0.36556
1999-04-01 0.74909
1999-07-01 0.74292
1999-10-01 0.73744
2000-01-01 0.98971
2000-04-01 0.78366
2000-07-01 0.91406
2000-10-01 0.71272
2001-01-01 0.95676
2001-04-01 0.70097
2001-07-01 0.28227
2001-10-01 -0.07487
2002-01-01 0.31944
2002-04-01 0.78622
2002-07-01 0.53826
2002-10-01 0.59136
2003-01-01 1.02865
2003-04-01 -0.16361
2003-07-01 0.74617
2003-10-01 0.37954
2004-01-01 0.84642
2004-04-01 0.78575
2004-07-01 0.63773
2004-10-01 1.07358
2005-01-01 0.50522
2005-04-01 0.67579
2005-07-01 1.51446
2005-10-01 0.93235
2006-01-01 0.52108
2006-04-01 0.90240
2006-07-01 0.94402
2006-10-01 -0.41050
2007-01-01 0.98056
2007-04-01 1.13255
2007-07-01 0.63301
2007-10-01 1.22680
2008-01-01 1.08319
2008-04-01 1.30094
2008-07-01 1.54172
2008-10-01 -2.29004
2009-01-01 -0.68787
2009-04-01 0.53160
2009-07-01 0.86039
2009-10-01 0.78293
2010-01-01 0.14330
2010-04-01 -0.00966
2010-07-01 0.30830
2010-10-01 0.76976
2011-01-01 1.03930
2011-04-01 1.21829
2011-07-01 0.65434
2011-10-01 0.38563
2012-01-01 0.52071
2012-04-01 0.34403
2012-07-01 0.42583
2012-10-01 0.59798
2013-01-01 0.29570
2013-04-01 0.10000
2013-07-01 0.53654
2013-10-01 0.28268
2014-01-01 0.47451
2014-04-01 0.74943
2014-07-01 0.27298
2014-10-01 -0.30127
DATE VALUE
1987-01-01 .
1987-04-01 3.32097
1987-07-01 1.91918
1987-10-01 1.31660
1988-01-01 0.74040
1988-04-01 3.89981
1988-07-01 1.77566
1988-10-01 1.02128
1989-01-01 1.69896
1989-04-01 2.71987
1989-07-01 1.10215
1989-10-01 0.19942
1990-01-01 0.27863
1990-04-01 1.11141
1990-07-01 -0.75896
1990-10-01 -1.64821
1991-01-01 -1.55517
1991-04-01 1.79763
1991-07-01 0.54849
1991-10-01 -0.67855
1992-01-01 -0.46885
1992-04-01 1.58816
1992-07-01 -0.10599
1992-10-01 -0.87533
1993-01-01 -0.37463
1993-04-01 1.36986
1993-07-01 0.76842
1993-10-01 -0.19721
1994-01-01 0.72454
1994-04-01 2.09260
1994-07-01 0.21778
1994-10-01 -0.43462
1995-01-01 -0.19258
1995-04-01 1.98096
1995-07-01 0.74420
1995-10-01 -0.45073
1996-01-01 0.12577
1996-04-01 1.88419
1996-07-01 0.75207
1996-10-01 -0.66079
1997-01-01 0.78837
1997-04-01 2.11440
1997-07-01 0.98145
1997-10-01 0.50966
1998-01-01 1.07311
1998-04-01 3.02182
1998-07-01 2.03851
1998-10-01 0.78801
1999-01-01 1.39852
1999-04-01 2.89965
1999-07-01 2.40633
1999-10-01 1.29857
2000-01-01 1.73975
2000-04-01 3.77000
2000-07-01 2.46699
2000-10-01 1.47654
2001-01-01 1.26969
2001-04-01 3.12986
2001-07-01 2.49357
2001-10-01 0.63203
2002-01-01 1.52284
2002-04-01 3.59322
2002-07-01 3.18226
2002-10-01 1.94244
2003-01-01 1.47768
2003-04-01 2.85101
2003-07-01 3.13711
2003-10-01 2.80327
2004-01-01 2.79008
2004-04-01 4.55353
2004-07-01 3.66858
2004-10-01 2.85750
2005-01-01 3.75935
2005-04-01 4.43880
2005-07-01 3.61064
2005-10-01 2.12475
2006-01-01 0.90389
2006-04-01 0.67317
2006-07-01 -0.95825
2006-10-01 -0.88778
2007-01-01 -0.86355
2007-04-01 -0.89812
2007-07-01 -1.72517
2007-10-01 -5.14416
2008-01-01 -6.67057
2008-04-01 -2.15236
2008-07-01 -3.49516
2008-10-01 -7.35646
2009-01-01 -7.34524
2009-04-01 3.11218
2009-07-01 3.24349
2009-10-01 -1.10537
2010-01-01 -2.87521
2010-04-01 4.69412
2010-07-01 -1.93086
2010-10-01 -3.48057
2011-01-01 -3.95752
2011-04-01 4.04900
2011-07-01 0.15291
2011-10-01 -3.80153
2012-01-01 -1.43628
2012-04-01 7.12503
2012-07-01 2.01413
2012-10-01 -0.44202
2013-01-01 1.09516
2013-04-01 7.20246
2013-07-01 3.09299
2013-10-01 -0.32452
2014-01-01 0.17276
## Prep workspace
options(width=120, digits=6, digits.secs=6)
rm(list=ls())
library(foreign)
library(grid)
library(plyr)
library(ggplot2)
library(tikzDevice)
print(options('tikzLatexPackages'))
options(tikzLatexPackages =
c("\\usepackage{tikz}\n",
"\\usepackage[active,tightpage,psfixbb]{preview}\n",
"\\PreviewEnvironment{pgfpicture}\n",
"\\setlength\\PreviewBorder{0pt}\n",
"\\usepackage{amsmath}\n",
"\\usepackage{xfrac}\n"
)
)
setTikzDefaults(overwrite = FALSE)
print(options('tikzLatexPackages'))
library(reshape)
library(vars)
library(scales)
library(zoo)
library(matpow)
## Define directories
scl.str.DAT_DIR <- "~/Dropbox/research/distant_speculators_and_mispricing/data/"
scl.str.FIG_DIR <- "~/Dropbox/research/distant_speculators_and_mispricing/figures/"
## Load CPI data
mat.df.CPI <- read.csv(file = paste(scl.str.DAT_DIR, "data--fred-quarterly-cpi.csv", sep = ""),
stringsAsFactors = FALSE
)
mat.df.CPI$DATE <- as.Date(mat.df.CPI$DATE, format = "%Y-%m-%d")
names(mat.df.CPI) <- c("t", "cpi")
scl.int.CPI_LEN <- dim(mat.df.CPI)[1]
mat.df.CPI$cpi <- as.numeric(mat.df.CPI$cpi)
mat.df.CPI$cpiLag <- c(NA, mat.df.CPI$cpi[1:(scl.int.CPI_LEN-1)])
mat.df.CPI <- mat.df.CPI[is.na(mat.df.CPI$cpiLag) == FALSE, ]
## Load HPI data
mat.df.HPI <- read.csv(file = paste(scl.str.DAT_DIR, "data--fred-quarterly-hpi.csv", sep = ""),
stringsAsFactors = FALSE
)
mat.df.HPI$DATE <- as.Date(mat.df.HPI$DATE, format = "%Y-%m-%d")
names(mat.df.HPI) <- c("t", "hpi")
scl.int.HPI_LEN <- dim(mat.df.HPI)[1]
mat.df.HPI$hpi <- as.numeric(mat.df.HPI$hpi)
mat.df.HPI$hpiLag <- c(NA, mat.df.HPI$hpi[1:(scl.int.HPI_LEN-1)])
mat.df.HPI <- mat.df.HPI[is.na(mat.df.HPI$hpiLag) == FALSE, ]
## Merge data
mat.df.DATA <- merge(mat.df.CPI, mat.df.HPI, by = c("t"))
## Demean variables
scl.flt.AVG_CPI <- mean(mat.df.DATA$cpi)
mat.df.DATA$cpi <- mat.df.DATA$cpi - scl.flt.AVG_CPI
mat.df.DATA$cpiLag <- mat.df.DATA$cpiLag - scl.flt.AVG_CPI
scl.flt.AVG_HPI <- mean(mat.df.DATA$hpi)
mat.df.DATA$hpi <- mat.df.DATA$hpi - scl.flt.AVG_HPI
mat.df.DATA$hpiLag <- mat.df.DATA$hpiLag - scl.flt.AVG_HPI
## Estimate CPI auto-regression
obj.lm.CPI_AR <- lm(cpi ~ 0 + cpiLag, data = mat.df.DATA)
summary(obj.lm.CPI_AR)
scl.flt.STD_EP <- sd(obj.lm.CPI_AR$residuals)
## Plot impulse response function for AR(1)
if (TRUE == TRUE) {
mat.df.PLOT <- data.frame(h = seq(0, 12),
dx = NA
)
for (h in 0:12) {
mat.df.PLOT$dx[h + 1] <- obj.lm.CPI_AR$coef^h * sd(obj.lm.CPI_AR$residuals)
}
theme_set(theme_bw())
scl.str.RAW_FILE <- "plot--cpi-ar1-irf"
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.AUX_FILE <- paste(scl.str.RAW_FILE,'.aux',sep='')
scl.str.LOG_FILE <- paste(scl.str.RAW_FILE,'.log',sep='')
tikz(file = scl.str.TEX_FILE, height = 2, width = 7, standAlone=TRUE)
obj.gg2.PLOT <- ggplot(data = mat.df.PLOT)
obj.gg2.PLOT <- obj.gg2.PLOT + scale_colour_brewer(palette="Set1")
obj.gg2.PLOT <- obj.gg2.PLOT + geom_hline(yintercept = 0,
linetype = 4,
size = 0.50
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = h,
y = dx
),
size = 1.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$h$ $\\scriptstyle [\\text{qtr}]$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\mathrm{Imp}(h)$ $\\scriptstyle [\\sfrac{\\%}{\\text{qtr}}]$")
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 12), breaks = c(0, 1, 2, 4, 8, 12))
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,0.15), "lines"),
axis.text = element_text(size = 8),
axis.title = element_text(size = 10),
plot.title = element_text(vjust = 1.75),
panel.grid.minor = element_blank()
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Impulse-Response Function for $\\mathrm{AR}(1)$")
print(obj.gg2.PLOT)
dev.off()
system(paste('lualatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('mv ', scl.str.PDF_FILE, ' ', scl.str.FIG_DIR, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
}
## Estimate CPI + HPI vector auto-regression
obj.lm.CPI_VAR <- lm(cpi ~ 0 + cpiLag + hpiLag, data = mat.df.DATA)
summary(obj.lm.CPI_VAR)
obj.lm.HPI_VAR <- lm(hpi ~ 0 + cpiLag + hpiLag, data = mat.df.DATA)
summary(obj.lm.HPI_VAR)
## Plot impulse response function for VAR
if (TRUE == TRUE) {
mat.flt.GAMMA <- matrix(rbind(obj.lm.CPI_VAR$coef, obj.lm.HPI_VAR$coef), nrow = 2)
mat.flt.SIGMA <- var(cbind(obj.lm.CPI_VAR$residuals, obj.lm.HPI_VAR$residuals))
mat.flt.C_TRN <- chol(mat.flt.SIGMA)
mat.flt.C <- t(chol(mat.flt.SIGMA))
## Correct
mat.df.PLOT1 <- data.frame(h = seq(0, 12),
dx = NA,
dy = NA
)
mat.df.PLOT1$dx[1] <- mat.flt.C[1,1]
mat.df.PLOT1$dy[1] <- mat.flt.C[2,1]
for (h in 1:12) {
mat.flt.GAMMA_POW <- matpow(mat.flt.GAMMA, h)$prod1
mat.df.PLOT1[h+1, c("dx", "dy")] <- mat.flt.GAMMA_POW %*% mat.flt.C %*% matrix(c(1,0), ncol = 1)
}
names(mat.df.PLOT1) <- c("h", "$\\Delta x_{t+h}$", "$\\Delta y_{t+h}$")
mat.df.PLOT1 <- melt(mat.df.PLOT1, c("h"))
names(mat.df.PLOT1) <- c("h", "variable", "vCorrect")
## Naive
mat.df.PLOT2 <- data.frame(h = seq(0, 12),
dx = NA,
dy = NA
)
mat.df.PLOT2$dx[1] <- sd(obj.lm.CPI_VAR$residuals)
mat.df.PLOT2$dy[1] <- 0
for (h in 1:12) {
mat.flt.GAMMA_POW <- matpow(mat.flt.GAMMA, h)$prod1
mat.df.PLOT2[h+1, c("dx", "dy")] <- mat.flt.GAMMA_POW %*% matrix(c(sd(obj.lm.CPI_VAR$residuals),0), ncol = 1)
}
names(mat.df.PLOT2) <- c("h", "$\\Delta x_{t+h}$", "$\\Delta y_{t+h}$")
mat.df.PLOT2 <- melt(mat.df.PLOT2, c("h"))
names(mat.df.PLOT2) <- c("h", "variable", "vNaive")
## Merge
mat.df.PLOT <- merge(mat.df.PLOT1, mat.df.PLOT2, by = c("h", "variable"))
mat.df.PLOT <- mat.df.PLOT[order(mat.df.PLOT$variable, mat.df.PLOT$h),]
theme_set(theme_bw())
scl.str.RAW_FILE <- "plot--cpi-vs-hpi-var-irf"
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.AUX_FILE <- paste(scl.str.RAW_FILE,'.aux',sep='')
scl.str.LOG_FILE <- paste(scl.str.RAW_FILE,'.log',sep='')
tikz(file = scl.str.TEX_FILE, height = 2, width = 7, standAlone=TRUE)
obj.gg2.PLOT <- ggplot(data = mat.df.PLOT)
obj.gg2.PLOT <- obj.gg2.PLOT + scale_colour_brewer(palette="Set1")
obj.gg2.PLOT <- obj.gg2.PLOT + geom_hline(yintercept = 0,
linetype = 4,
size = 0.50
)
## obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = h,
## y = vNaive,
## group = variable
## ),
## size = 1.25,
## colour = "red"
## )
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = h,
y = vCorrect,
group = variable
),
size = 1.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$h$ $\\scriptstyle [\\text{qtr}]$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\mathrm{Imp}_x(h)$ $\\scriptstyle [\\sfrac{\\%}{\\text{qtr}}]$")
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_continuous(limits = c(0, 12), breaks = c(0, 1, 2, 4, 8, 12))
obj.gg2.PLOT <- obj.gg2.PLOT + facet_wrap(~ variable, nrow = 1)
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,0.15,0.15,0.15), "lines"),
axis.text = element_text(size = 8),
axis.title = element_text(size = 10),
plot.title = element_text(vjust = 1.75),
panel.grid.minor = element_blank()
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Impulse-Response Function for $\\mathrm{VAR}$")
print(obj.gg2.PLOT)
dev.off()
system(paste('lualatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('mv ', scl.str.PDF_FILE, ' ', scl.str.FIG_DIR, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment