Skip to content

Instantly share code, notes, and snippets.

@alexchinco

alexchinco/asset-pricing-time-scales.R

Created August 22, 2013 16:41
Show Gist options
  • Save alexchinco/6309722 to your computer and use it in GitHub Desktop.
Save alexchinco/6309722 to your computer and use it in GitHub Desktop.
Download ZIP
Code to create figures for my post "Identifying Relevant Asset Pricing Time Scales". url: http://www.alexchinco.com/identifying-relevant-asset-pricing-time-scales/
Raw
asset-pricing-time-scales.R
################################################################################
################################################################################
## Prep workspace
################################################################################
################################################################################
options(width=200, digits=6, digits.secs=6)
rm(list=ls())
library(foreign)
library(grid)
library(plyr)
library(ggplot2)
library(tikzDevice)
library(reshape)
library(vars)
library(scales)
library(dynlm)
scl.str.DAT_DIR <- "~/Dropbox/research/technical_notes/"
scl.str.FIG_DIR <- "~/Dropbox/research/technical_notes/"
set.seed(1234)
################################################################################
################################################################################
## Create figures for sections 2, 3, and 4
################################################################################
################################################################################
################################################################################
## Simulation log price time series
################################################################################
scl.int.T <- 365 * 20
vec.flt.X <- rep(NA, scl.int.T)
vec.flt.dX <- rep(NA, scl.int.T)
scl.int.PERIOD_1 <- 7
scl.int.PERIOD_2 <- 30
scl.flt.THETA_1 <- 95/100
scl.flt.THETA_2 <- -95/100
scl.flt.SIGMA <- 1/10
vec.flt.X[1:100] <- scl.flt.SIGMA * rnorm(100,0,1)
for (t in 101:scl.int.T) {
vec.flt.X[t] <- scl.flt.THETA_1 * mean(vec.flt.X[(t-scl.int.PERIOD_1):(t-1)]) + scl.flt.THETA_2 * mean(vec.flt.X[(t-scl.int.PERIOD_2):(t-1)]) + scl.flt.SIGMA * rnorm(1,0,1)
}
mat.dfm.DATA <- data.frame(t = seq(1, scl.int.T),
tD = seq(as.Date("1980-01-01"), len = scl.int.T, by = "1 day"),
dx = c(NA, vec.flt.X[-1] - vec.flt.X[-scl.int.T]),
x = vec.flt.X
)
################################################################################
## Plot return time series
################################################################################
theme_set(theme_bw())
mat.dfm.PLOT <- mat.dfm.DATA[mat.dfm.DATA$t %in% seq(10 * 365 + 1, 11 * 365), c("tD", "dx", "x")]
names(mat.dfm.PLOT) <- c("tD", "$r_t$", "$\\log p_t$")
mat.dfm.PLOT <- melt(mat.dfm.PLOT, c("tD"))
scl.str.RAW_FILE <- 'plot-multiscale-process'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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.dfm.PLOT)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = tD,
y = value,
group = variable,
colour = variable
),
size = 0.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$t \\ {\\scriptstyle (\\mathrm{days})}$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
obj.gg2.PLOT <- obj.gg2.PLOT + coord_cartesian(ylim = c(-0.60, 0.60))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(breaks = c(-0.50, 0, 0.50))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_date(labels = date_format("%b"), breaks = date_breaks("months"))
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = as.Date("1990-09-01"),
y = -0.50,
size = 3,
label = "$\\log p_t = \\frac{95}{100} \\cdot \\left( \\frac{1}{7} \\cdot \\sum_{h=1}^h \\log p_{t - h{\\scriptscriptstyle \\mathrm{days}}} \\right) - \\frac{95}{100} \\cdot \\left( \\frac{1}{30} \\cdot \\sum_{h=1}^{30} \\log p_{t - h{\\scriptscriptstyle \\mathrm{days}}} \\right) + \\frac{1}{10} \\cdot \\varepsilon_t$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Multi-Scale Log Price Process")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,1,0.0,0.0), "lines"),
plot.title = element_text(vjust = 1.75),
legend.position = c(0.05,0.85),
legend.title = element_text(colour = NA, size = 0),
legend.text = element_text(size = 8),
legend.background = element_rect(fill=alpha('blue', 0))
)
print(obj.gg2.PLOT)
dev.off()
system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
################################################################################
## Run predictive regression
################################################################################
mat.dfm.DATA$L0 <- mat.dfm.DATA$x
scl.str.FORMULA <- 'x ~ 1'
vec.int.LAGS <- seq(1, len = 90, by = 1)
for (h in vec.int.LAGS) {
mat.dfm.DATA[, paste('L', h, sep = '')] <- c(NA, mat.dfm.DATA[1:(scl.int.T - 1), paste('L', h-1, sep = '')])
}
mat.dfm.AC <- data.frame(n = vec.int.LAGS,
h = vec.int.LAGS,
c = NA,
ub = NA,
lb = NA
)
for (n in vec.int.LAGS) {
obj.lm.AC <- lm(paste(scl.str.FORMULA, ' + L', n, sep = ''), data = mat.dfm.DATA[91:scl.int.T,])
mat.dfm.AC[mat.dfm.AC$n == n, ]$c <- summary(obj.lm.AC)$coef[2,1]
mat.dfm.AC[mat.dfm.AC$n == n, ]$ub <- summary(obj.lm.AC)$coef[2,1] + 2.64 * summary(obj.lm.AC)$coef[2,2]
mat.dfm.AC[mat.dfm.AC$n == n, ]$lb <- summary(obj.lm.AC)$coef[2,1] - 2.64 * summary(obj.lm.AC)$coef[2,2]
}
theme_set(theme_bw())
scl.str.RAW_FILE <- 'plot-multiscale-autoregression-coefficients'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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.dfm.AC)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_hline(yintercept = 0,
colour = "black",
size = 0.25,
linetype = 4
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = h,
y = c
),
colour = "black",
size = 0.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$h \\ {\\scriptstyle (\\mathrm{days})}$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\hat{\\mathrm{C}}(h)$")
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 75,
y = 0.50,
size = 3,
label = "$\\log p_t = \\hat{\\mathrm{C}}(h) \\cdot \\log p_{t-h} + \\xi_{h,t}$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Autocorrelation Function")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,1,0.0,0.0), "lines"),
plot.title = element_text(vjust = 1.75)
)
print(obj.gg2.PLOT)
dev.off()
system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
################################################################################
## Run frequency regression
################################################################################
vec.flt.FREQS <- seq(1/100, 1/3, length.out = 250)
vec.flt.TIME <- seq(1, scl.int.T)
scl.int.N <- length(vec.flt.FREQS)
mat.dfm.PS <- data.frame(n = 1:scl.int.N,
f = vec.flt.FREQS,
a = NA,
b = NA
)
for (n in 1:scl.int.N) {
obj.lm.PS <- lm(mat.dfm.DATA$x ~ 0 + sin(vec.flt.FREQS[n] * vec.flt.TIME) + cos(vec.flt.FREQS[n] * vec.flt.TIME))
mat.dfm.PS[mat.dfm.PS$n == n, ]$a <- summary(obj.lm.PS)$coef[1,1]
mat.dfm.PS[mat.dfm.PS$n == n, ]$b <- summary(obj.lm.PS)$coef[2,1]
}
theme_set(theme_bw())
mat.dfm.PS$s <- (mat.dfm.PS$a^2 + mat.dfm.PS$b^2)/2
scl.str.RAW_FILE <- 'plot-multiscale-process-powerspectrum'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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.dfm.PS)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_hline(yintercept = 0,
colour = "black",
size = 0.25,
linetype = 4
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = f,
y = log10(s)
),
colour = "black",
size = 0.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.25,
y = -1,
size = 3,
label = "$\\log p_t = \\hat{a}_f \\cdot \\sin(f \\cdot t) + \\hat{b}_f \\cdot \\cos( f \\cdot t) + \\xi_{f,t}$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.25,
y = -2.5,
size = 3,
label = "$\\hat{\\mathrm{S}}(f) = \\frac{1}{2} \\cdot (\\hat{a}_f^2 + \\hat{b}_f^2)$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(labels = math_format(10^.x)) + annotation_logticks(sides = "l")
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$f \\ {\\scriptstyle (1/\\mathrm{days})}$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\hat{\\mathrm{S}}(f)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Power Spectrum")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,1,0.0,0.0), "lines"),
plot.title = element_text(vjust = 1.75)
)
print(obj.gg2.PLOT)
dev.off()
system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
################################################################################
## Recover autocorrelation coefficients from power spectrum
################################################################################
mat.dfm.AC$cHat <- NA
for (n in 1:dim(mat.dfm.AC)[1]) {
mat.dfm.AC$cHat[n] <- sum(mat.dfm.PS$s * cos(mat.dfm.PS$f * mat.dfm.AC$h[n]))/sum(mat.dfm.PS$s)
}
obj.lm.PvA <- lm(cHat ~ c, data = mat.dfm.AC)
theme_set(theme_bw())
scl.str.RAW_FILE <- 'plot-multiscale-actual-vs-predicted-autoregression-coefficients'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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.dfm.AC)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_abline(yintercept = obj.lm.PvA$coef[1],
slope = obj.lm.PvA$coef[2],
colour = "red",
size = 0.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_point(aes(x = c,
y = cHat
),
colour = "black",
size = 0.50
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = -0.50,
y = 0.75,
size = 3,
label = "$\\hat{\\mathrm{C}}_{WK}(h) = \\sum_f \\left( \\frac{\\hat{\\mathrm{S}}(f)}{\\sum_{f'} \\hat{\\mathrm{S}}(f')} \\right) \\cdot \\cos(f \\cdot h)$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.725,
y = 0.25,
size = 3,
label = "$\\hat{\\mathrm{C}}_{WK}(h) = \\alpha + \\beta \\cdot \\hat{\\mathrm{C}}(h) + \\epsilon_h$",
colour = "red"
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.655,
y = 0.00,
size = 3,
label = paste("$\\alpha = ", round(obj.lm.PvA$coef[1],2),"$", sep = ""),
colour = "red"
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.675,
y = -0.25,
size = 3,
label = paste("$\\beta = ", round(obj.lm.PvA$coef[2],2),"$", sep = ""),
colour = "red"
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.685,
y = -0.50,
size = 3,
label = paste("$R^2 = ", round(summary(obj.lm.PvA)$r.squared * 100,3),"{\\scriptstyle \\%}$", sep = ""),
colour = "red"
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$\\hat{\\mathrm{C}}(h)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\hat{\\mathrm{C}}_{WK}(h)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Empirical vs Wiener-Khintchine Implied Autocorrelation Function")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,1,0.0,0.0), "lines"),
plot.title = element_text(vjust = 1.75)
)
print(obj.gg2.PLOT)
dev.off()
system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
################################################################################
################################################################################
## Create figures for section 5
################################################################################
################################################################################
################################################################################
## Simulation log price time series
################################################################################
scl.int.T <- 365 * 20
vec.flt.X_1 <- rep(NA, scl.int.T)
vec.flt.dX_1 <- rep(NA, scl.int.T)
vec.flt.X_2 <- rep(NA, scl.int.T)
vec.flt.dX_2 <- rep(NA, scl.int.T)
vec.flt.X_3 <- rep(NA, scl.int.T)
vec.flt.dX_3 <- rep(NA, scl.int.T)
vec.flt.X_4 <- rep(NA, scl.int.T)
vec.flt.dX_4 <- rep(NA, scl.int.T)
scl.int.PERIOD_1 <- 30
scl.int.PERIOD_2 <- 7
scl.int.PERIOD_3 <- 1
scl.int.PERIOD_4 <- 1/24
scl.flt.THETA_1 <- 1 - exp(-log(2)/(scl.int.PERIOD_1/2))
scl.flt.THETA_2 <- 1 - exp(-log(2)/(scl.int.PERIOD_2/2))
scl.flt.THETA_3 <- 1 - exp(-log(2)/(scl.int.PERIOD_3/2))
scl.flt.THETA_4 <- 1 - exp(-log(2)/(scl.int.PERIOD_4/2))
scl.flt.SIGMA <- 1/10
vec.flt.X_1[1] <- vec.flt.X_2[1] <- vec.flt.X_3[1] <- vec.flt.X_4[1] <- 0
for (t in 2:scl.int.T) {
vec.flt.dX_1[t] <- - scl.flt.THETA_1 * vec.flt.X_1[t-1] + scl.flt.SIGMA * rnorm(1,0,1)
vec.flt.X_1[t] <- vec.flt.X_1[t-1] + vec.flt.dX_1[t]
vec.flt.dX_2[t] <- - scl.flt.THETA_2 * vec.flt.X_2[t-1] + scl.flt.SIGMA * rnorm(1,0,1)
vec.flt.X_2[t] <- vec.flt.X_2[t-1] + vec.flt.dX_2[t]
vec.flt.dX_3[t] <- - scl.flt.THETA_3 * vec.flt.X_3[t-1] + scl.flt.SIGMA * rnorm(1,0,1)
vec.flt.X_3[t] <- vec.flt.X_3[t-1] + vec.flt.dX_3[t]
vec.flt.dX_4[t] <- - scl.flt.THETA_4 * vec.flt.X_4[t-1] + scl.flt.SIGMA * rnorm(1,0,1)
vec.flt.X_4[t] <- vec.flt.X_4[t-1] + vec.flt.dX_4[t]
}
mat.dfm.DATA <- data.frame(t = seq(1, scl.int.T),
tD = seq(as.Date("1980-01-01"), len = scl.int.T, by = "1 day"),
dx1 = vec.flt.dX_1,
x1 = vec.flt.X_1,
dx2 = vec.flt.dX_2,
x2 = vec.flt.X_2,
dx3 = vec.flt.dX_3,
x3 = vec.flt.X_3,
dx4 = vec.flt.dX_4,
x4 = vec.flt.X_4
)
################################################################################
## Plot return time series
################################################################################
theme_set(theme_bw())
mat.dfm.PLOT <- mat.dfm.DATA[mat.dfm.DATA$t %in% seq(10 * 365 + 1, 11 * 365), c("tD", "dx1", "x1")]
names(mat.dfm.PLOT) <- c("tD", "$r_t$", "$\\log p_t$")
mat.dfm.PLOT <- melt(mat.dfm.PLOT, c("tD"))
scl.str.RAW_FILE <- 'plot-ou-process'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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.dfm.PLOT)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = tD,
y = value,
group = variable,
colour = variable
),
size = 0.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$t \\ {\\scriptstyle (\\mathrm{days})}$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("")
obj.gg2.PLOT <- obj.gg2.PLOT + coord_cartesian(ylim = c(-0.6, 0.6))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(breaks = c(-0.5, 0, 0.5))
obj.gg2.PLOT <- obj.gg2.PLOT + scale_x_date(labels = date_format("%b"), breaks = date_breaks("months"))
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = as.Date("1990-04-15"),
y = 0.50,
size = 3,
label = paste("$\\log p_{t + 1{\\scriptscriptstyle \\mathrm{day}}} = \\log p_t - \\left(1 - e^{-\\log(2)/", round(scl.int.PERIOD_1/2,1),"} \\right) \\cdot \\log p_t \\cdot 1{\\scriptstyle \\mathrm{day}} + \\frac{1}{10} \\cdot \\varepsilon_{t+1{\\scriptscriptstyle \\mathrm{day}}}$", sep = "")
)
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Ornstein-Uhlenbeck Process")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,1,0.0,0.0), "lines"),
plot.title = element_text(vjust = 1.75),
legend.position = c(0.10,0.20),
legend.title = element_text(colour = NA, size = 0),
legend.text = element_text(size = 8),
legend.background = element_rect(fill=alpha('blue', 0))
)
print(obj.gg2.PLOT)
dev.off()
system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, sep = ''))
system(paste('rm ', scl.str.AUX_FILE, sep = ''))
system(paste('rm ', scl.str.LOG_FILE, sep = ''))
################################################################################
## Run frequency regression
################################################################################
vec.flt.FREQS <- seq(1/250, 1, length.out = 250)
vec.flt.TIME <- seq(1, scl.int.T)
scl.int.N <- length(vec.flt.FREQS)
mat.dfm.PS <- data.frame(n = 1:scl.int.N,
f = vec.flt.FREQS,
a1 = NA,
b1 = NA,
a2 = NA,
b2 = NA,
a3 = NA,
b3 = NA,
a4 = NA,
b4 = NA
)
for (n in 1:scl.int.N) {
obj.lm.PS <- lm(mat.dfm.DATA$x1 ~ 0 + sin(vec.flt.FREQS[n] * vec.flt.TIME) + cos(vec.flt.FREQS[n] * vec.flt.TIME))
mat.dfm.PS[mat.dfm.PS$n == n, ]$a1 <- summary(obj.lm.PS)$coef[1,1]
mat.dfm.PS[mat.dfm.PS$n == n, ]$b1 <- summary(obj.lm.PS)$coef[2,1]
obj.lm.PS <- lm(mat.dfm.DATA$x2 ~ 0 + sin(vec.flt.FREQS[n] * vec.flt.TIME) + cos(vec.flt.FREQS[n] * vec.flt.TIME))
mat.dfm.PS[mat.dfm.PS$n == n, ]$a2 <- summary(obj.lm.PS)$coef[1,1]
mat.dfm.PS[mat.dfm.PS$n == n, ]$b2 <- summary(obj.lm.PS)$coef[2,1]
obj.lm.PS <- lm(mat.dfm.DATA$x3 ~ 0 + sin(vec.flt.FREQS[n] * vec.flt.TIME) + cos(vec.flt.FREQS[n] * vec.flt.TIME))
mat.dfm.PS[mat.dfm.PS$n == n, ]$a3 <- summary(obj.lm.PS)$coef[1,1]
mat.dfm.PS[mat.dfm.PS$n == n, ]$b3 <- summary(obj.lm.PS)$coef[2,1]
obj.lm.PS <- lm(mat.dfm.DATA$x4 ~ 0 + sin(vec.flt.FREQS[n] * vec.flt.TIME) + cos(vec.flt.FREQS[n] * vec.flt.TIME))
mat.dfm.PS[mat.dfm.PS$n == n, ]$a4 <- summary(obj.lm.PS)$coef[1,1]
mat.dfm.PS[mat.dfm.PS$n == n, ]$b4 <- summary(obj.lm.PS)$coef[2,1]
}
theme_set(theme_bw())
mat.dfm.PS$s1 <- (mat.dfm.PS$a1^2 + mat.dfm.PS$b1^2)/2
mat.dfm.PS$s2 <- (mat.dfm.PS$a2^2 + mat.dfm.PS$b2^2)/2
mat.dfm.PS$s3 <- (mat.dfm.PS$a3^2 + mat.dfm.PS$b3^2)/2
mat.dfm.PS$s4 <- (mat.dfm.PS$a4^2 + mat.dfm.PS$b4^2)/2
mat.dfm.PLOT <- mat.dfm.PS[, c("f", "s1", "s2", "s3", "s4")]
names(mat.dfm.PLOT) <- c("f", "Month", "Week", "Day", "Hour")
mat.dfm.PLOT <- melt(mat.dfm.PLOT, c("f"))
names(mat.dfm.PLOT) <- c("f", "Period", "value")
scl.str.RAW_FILE <- 'plot-ou-process-powerspectrum'
scl.str.TEX_FILE <- paste(scl.str.RAW_FILE,'.tex',sep='')
scl.str.PDF_FILE <- paste(scl.str.RAW_FILE,'.pdf',sep='')
scl.str.PNG_FILE <- paste(scl.str.RAW_FILE,'.png',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.dfm.PLOT)
obj.gg2.PLOT <- obj.gg2.PLOT + geom_path(aes(x = f,
y = log10(value),
group = Period,
colour = Period
),
size = 0.25
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.75,
y = -3,
size = 3,
label = "$\\log p_t = \\hat{a}_f \\cdot \\sin(f \\cdot t) + \\hat{b}_f \\cdot \\cos( f \\cdot t) + \\xi_{f,t}$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + annotate("text",
x = 0.75,
y = -4,
size = 3,
label = "$\\hat{\\mathrm{S}}(f) = \\frac{1}{2} \\cdot (\\hat{a}_f^2 + \\hat{b}_f^2)$"
)
obj.gg2.PLOT <- obj.gg2.PLOT + scale_y_continuous(labels = math_format(10^.x)) + annotation_logticks(sides = "l")
obj.gg2.PLOT <- obj.gg2.PLOT + xlab("$f \\ {\\scriptstyle (1/\\mathrm{days})}$")
obj.gg2.PLOT <- obj.gg2.PLOT + ylab("$\\hat{\\mathrm{S}}(f)$")
obj.gg2.PLOT <- obj.gg2.PLOT + ggtitle("Power Spectrum")
obj.gg2.PLOT <- obj.gg2.PLOT + theme(plot.margin = unit(c(1,1,0.0,0.0), "lines"),
plot.title = element_text(vjust = 1.75),
legend.direction = "horizontal",
legend.position = c(0.25,0.15),
legend.title = element_text(size = 8),
legend.text = element_text(size = 8)
)
print(obj.gg2.PLOT)
dev.off()
system(paste('pdflatex', file.path(scl.str.TEX_FILE)), ignore.stdout = TRUE)
system(paste('rm ', scl.str.TEX_FILE, 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