Nothing
##### TIMSAC84 #####
baysea <- function(y, period = 12, span = 4, shift = 1, forecast = 0,
trend.order = 2, seasonal.order = 1, year = 0, month = 1,
out = 0, rigid = 1, zersum = 1, delta = 7, alpha = 0.01,
beta = 0.01, gamma = 0.1, spec = TRUE, plot = TRUE,
separate.graphics = FALSE)
{
if (seasonal.order > span)
stop("seasonal.order is smaller than or equal to span")
if (span < 1)
stop("span is greater than or equal to 1")
if (trend.order < 1)
stop("trend.order is greater than or equal to 1")
ndata <- length(y)
npf <- ndata + forecast
ipara <- rep(0, 12)
ipara[1] <- period
ipara[2] <- span
ipara[3] <- shift
ipara[4] <- trend.order
ipara[5] <- seasonal.order
ipara[6] <- 0 # logt
ipara[7] <- year
ipara[8] <- month
ipara[9] <- 1 # nday
if (spec == TRUE)
ipara[10] <- 1 # spectrum estimation option
if (spec == FALSE)
ipara[10] <- 0
ipara[11] <- out # ioutd : outlier correction option
para <- rep(0, 8)
para[1] <- rigid # controls the rigidity of the seasonal component
para[2] <- 1 # wtrd
para[3] <- 1 # dd
para[4] <- zersum # controls the sum of the seasonals within a period
para[5] <- delta # controls the leap year effect
para[6] <- alpha # controls prior variance of initial trend
para[7] <- beta # controls prior variance of initial seasonal
para[8] <- gamma # controls prior variance of initial sum of seasonal
arft <- rep(0, 3)
arfs <- rep(0, 3)
arfn <- rep(0, 3)
# subroutine arcoef ---> subroutine partar
iart <- 0
# for(i in 1:3) if(arft[i] != 0.0) iord <- i
# if(iord != 0) {
# arm <- matrix(0, dim = c(3,3))
# for(i in 1:3) arm[i,i] <- arft[i]
# arm[2,1] <- arm[1,1]-arft[2]*arm[1,1]
# arm[3,1] <- arm[2,1]-arft[3]*arm[2,2]
# arm[3,2] <- arm[2,2]-arft[3]*arm[2,1]
# for(i in 1:iord) arft[i] <- arm[iord,i]
# }
iars <- 0
iarn <- 0
is <- period * seasonal.order
iprd <- 2
if (period == 1)
iprd <- 1
lftrn <- trend.order + iart
lfsea <- (seasonal.order + iars) * period + iarn
idc <- lftrn * iprd + 1
idcx <- lfsea * 2 + 1
if ((period > 1) && (idc < idcx))
idc <- idcx
if ((period > 1) && (idc < period*2-1))
idc <- period * 2 - 1
ipara[12] <- idc
z <- .Fortran(C_bayseaf,
as.double(y),
as.integer(ndata),
as.integer(forecast),
outlier = double(ndata),
dmoi = double(ndata),
trend = double(npf),
season = double(npf),
tdcmp = double(npf),
irreg = double(ndata),
adjust = double(ndata),
est = double(npf),
psds = double(npf),
psdt = double(npf),
avabic = double(1),
as.integer(ipara),
as.double(para),
as.double(arft),
as.double(arfs),
as.double(arfn),
as.integer(iart),
as.integer(iars),
as.integer(iarn))
outlier <- NULL
if (out != 0)
outlier <- z$outlier
tre <- z$trend
sea <- z$season
tday <- NULL
if (year != 0)
tday <- z$tdcmp
irr <- z$irreg
adj <- z$adjust
est <- z$est
psds <- z$psds
psdt <- z$psdt
abic <- z$avabic
baysea.out <- list(outlier = outlier, trend = tre, season = sea, tday = tday,
irregular = irr, adjust = adj, smoothed = est,
aveABIC = abic)
if (spec == TRUE) {
z1 <- bayspec(y, period, tre, trend.order, sea, seasonal.order, irr, adj)
spec1 <- z1$irregular.spec
spec2 <- z1$adjusted.spec
spec3 <- z1$differenced.trend
spec4 <- z1$differenced.season
if (plot == TRUE)
plotCompSpec(y, period, outlier, tre, trend.order, sea,
seasonal.order, tday, irr, adj, est, psdt, psds, spec1,
spec2, spec3, spec4, spec, separate.graphics)
ir.spec <- list(acov = spec1$acov, acor = spec1$acor, mean = spec1$mean,
v = spec1$v, aic = spec1$aic, parcor = spec1$parcor,
rspec = spec1$rspec)
ad.spec <- list(acov = spec2$acov, acor = spec2$acor, mean = spec2$mean,
v = spec2$v, aic = spec2$aic, parcor = spec2$parcor,
rspec = spec2$rspec)
diff.trend <- list(acov = spec3$acov, acor = spec3$acor,
mean = spec3$mean, v = spec3$v, aic = spec3$aic,
parcor = spec3$parcor)
diff.season <- list(acov = spec4$acov, acor = spec4$acor,
mean = spec4$mean, v = spec4$v, aic = spec4$aic,
parcor = spec4$parcor)
baysea.out <- c(baysea.out, irregular.spec = ir.spec,
adjusted.spec = ad.spec, differenced.trend = diff.trend,
differenced.season = diff.season)
} else {
if (plot == TRUE) {
ir.spec <- NULL
ad.spec <- NULL
diff.trend <- NULL
diff.season <- NULL
plotCompSpec(y, period, outlier, tre, trend.order, sea,
seasonal.order, tday, irr, adj, est, psdt, psds, ir.spec,
ad.spec, diff.trend, diff.season, spec,
separate.graphics)
}
}
return(baysea.out)
}
bayspec <- function(y, period, trend, trend.order, season, seasonal.order,
irregular, adjust)
{
ndata <- length(y)
npf <- length(trend)
lag <- min((ndata-1), 60)
lag1 <- lag + 1
# ifpl <- min(3*sqrt(ndata), 50, lag)
ifpl <- min(30, ndata-1)
ifpl1 <- ifpl + 1
# SPECTRUM OF IRREGULAR
mode <- 1
z <- .Fortran(C_spgrh,
as.double(irregular),
as.integer(ndata),
as.integer(lag1),
as.integer(ifpl1),
as.integer(mode),
as.integer(period),
cxx = double(lag1),
cn = double(lag1),
xmean = double(1),
sd = double(ifpl1),
aic = double(ifpl1),
parcor = double(ifpl),
pxx = double(lag1),
ier = integer(1))
pxx <- z$pxx
sxx <- rep(0, lag1)
for(i in 1:lag1) {
t <- abs(pxx[i])
sxx[i] <- log10(t)
}
if (z$ier == 2600)
warning("accuracy of computation lost")
irregular.spec <- list(n = ndata, acov = z$cxx, acor = z$cn,
mean = z$xmean, v = z$sd, aic = z$aic,
parcor = z$parcor, rspec = pxx, rpspec = sxx)
# SPECTRUM OF DIFFERENCED ADJUSTED SERIES
n1 <- ndata - 1
dadj <- adjust
for(i in 1:n1)
dadj[i] <- dadj[i+1] - dadj[i]
dadj <- dadj[1:n1]
lag <- min((n1-1), 60)
lag1 <- lag + 1
z <- .Fortran(C_spgrh,
as.double(dadj),
as.integer(n1),
as.integer(lag1),
as.integer(ifpl1),
as.integer(mode),
as.integer(period),
cxx = double(lag1),
cn = double(lag1),
xmean = double(1),
sd = double(ifpl1),
aic = double(ifpl1),
parcor = double(ifpl),
pxx = double(lag1),
ier = integer(1))
pxx <- z$pxx
sxx <- rep(0, lag1)
for(i in 1:lag1) {
t <- abs(pxx[i])
sxx[i] <- log10(t)
}
if (z$ier == 2600)
cat(" ***** Warnning : Accuracy of computation lost\n")
adjusted.spec <- list(n = n1, acov = z$cxx, acor = z$cn, mean = z$xmean,
v = z$sd, aic = z$aic, parcor = z$parcor,
rspec = pxx, rpspec = sxx)
# PARCOR OF TREND.ORDER TIME(S) DIFFERENCED TREND SERIES
n1 <- ndata
trendd <- trend
for(j in 1:trend.order) {
n1 <- n1 - 1
for(i in 1:n1)
trendd[i] <- trendd[i+1] - trendd[i]
}
trendd <- trendd[1:n1]
lag <- min((n1-1), 60)
lag1 <- lag + 1
mode <- 0
z <- .Fortran(C_spgrh,
as.double(trendd),
as.integer(n1),
as.integer(lag1),
as.integer(ifpl1),
as.integer(mode),
as.integer(period),
cxx = double(lag1),
cn = double(lag1),
xmean = double(1),
sd = double(ifpl1),
aic = double(ifpl1),
parcor = double(ifpl),
pxx = double(lag1),
ier = integer(1))
if (z$ier == 2600)
warning("accuracy of computation lost")
differenced.trend <- list(n = n1, acov = z$cxx, acor = z$cn,
mean = z$xmean, v = z$sd, aic = z$aic,
parcor = z$parcor)
# PARCOR OF SEASONAL.ORDER TIME(S) DIFFERENCED SEASONAL SERIES
n1 <- ndata
seasond <- season
for(j in 1:seasonal.order) {
n1 <- n1 - period
for(i in 1:n1)
seasond[i] <- seasond[i+period] - seasond[i]
}
seasond <- seasond[1:n1]
lag <- min((n1-1), 60)
lag1 <- lag + 1
z <- .Fortran(C_spgrh,
as.double(seasond),
as.integer(n1),
as.integer(lag1),
as.integer(ifpl1),
as.integer(mode),
as.integer(period),
cxx = double(lag1),
cn = double(lag1),
xmean = double(1),
sd = double(ifpl1),
aic = double(ifpl1),
parcor = double(ifpl),
pxx = double(lag1),
ier = integer(1))
if (z$ier == 2600)
warning("accuracy of computation lost")
differenced.season <- list(n = n1, acov = z$cxx, acor = z$cn,
mean = z$xmean, v = z$sd, aic = z$aic,
parcor = z$parcor)
bayspec.out <- list(irregular.spec = irregular.spec,
adjusted.spec = adjusted.spec,
differenced.trend = differenced.trend,
differenced.season = differenced.season)
}
plotCompSpec <- function(y, period, outlier, trend, trend.order, season,
seasonal.order, tday, irregular, adjust, smoothed, psdt,
psds, irregular.spec, adjusted.spec, differenced.trend,
differenced.season, spec, separate.graphics)
{
oldpar <- par(no.readonly = TRUE)
ndata <- length(y)
npf <- length(trend)
ymax1 <- max(y, trend, adjust, smoothed)
ymin1 <- min(y, trend, adjust, smoothed)
plot(y, type = "l", main = "Original Data", xlab = "time", xlim = c(0,npf),
ylim = c(ymin1,ymax1))
nw <- 1
if (separate.graphics == TRUE) {
dev.new()
} else {
par(ask = TRUE)
}
plot(trend, type = "l", main = "Trend and 2*(post SD)", cex.main = 0.9,
xlab = "time", ylab = "", xlim = c(0, npf), ylim = c(ymin1, ymax1))
par(new = TRUE)
xtem <- trend + psdt
plot(xtem, type = "l", lty = 3, main = "", xlab = "", ylab = "",
xlim = c(0, npf), ylim = c(ymin1, ymax1))
par(new = TRUE)
xtem <- trend - psdt
plot(xtem, type = "l", lty = 3, main = "", xlab = "", ylab = "",
xlim = c(0, npf), ylim = c(ymin1, ymax1))
if (separate.graphics == TRUE)
dev.new()
plot(adjust, pch = 0, type = "l",
main = "Adjusted = Original Data - Seasonal - Trading.Day.Comp - Outlier",
cex.main = 0.9, xlab = "time", ylab = "", xlim = c(0, npf),
ylim = c(ymin1, ymax1))
if (separate.graphics == TRUE)
dev.new()
plot(smoothed, pch = 0, type = "l",
main = "Smoothed = Trend + Seasonal + Trading.Day.Comp", cex.main = 0.9,
xlab = "time", ylab = "", xlim = c(0, npf), ylim = c(ymin1, ymax1))
ymax2 <- max(irregular)
ymin2 <- min(irregular)
if (seasonal.order != 0) {
ymax2 <- max(season, ymax2)
ymin2 <- min(season, ymin2)
}
if (is.null(tday) == FALSE) {
ymax2 <- max(tday, ymax2)
ymin2 <- min(tday, ymin2)
}
my <- max(ymax2, abs(ymin2)) * 1.5
if (seasonal.order != 0) {
if (separate.graphics == TRUE)
dev.new()
plot(season, type = "l", main = "Seasonal and 2*(post SD)",
cex.main = 0.9, xlab = "time", ylab = "", ylim = c(-my, my))
par(new = TRUE)
xtem <- season + psds
plot(xtem, type = "l", lty = 3, main = "", xlab = "", ylab = "",
ylim = c(-my, my))
par(new = TRUE)
xtem <- season - psds
plot(xtem, type = "l", lty = 3, main = "", xlab = "", ylab = "",
ylim = c(-my, my))
}
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(2, 1))
plot(irregular, type = "l",
main = "Irregular = Original Data - Trend - Seasonal - Trading.Day.Comp",
cex.main = 0.9, xlab = "time", ylab = "", ylim = c(-my, my))
vy <- sd(irregular) * 5
plot(irregular, type = "l",
main = "Irregular ( Scaled by the Standard Deviation)", cex.main = 0.9,
xlab = "time", ylab = "", ylim = c(-vy, vy))
par(mfrow = c(1, 1))
# SPECTRUM OF IRREGULAR
if (spec == TRUE) {
lag <- length(irregular.spec$acov) - 1
n <- irregular.spec$n
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(3, 1))
plot((0:lag), irregular.spec$acor, type = "h",
main = "Autocorrelation & Parcor of Irregular",
ylab = "Autocorrelation", xlab = "lag", ylim = c(-1, 1))
plot(irregular.spec$parcor, type = "h", ylab = "Parcor", xlab = "lag",
ylim = c(-1, 1))
abline(h = 0, lty = 1)
abline(h = 2/sqrt(n), lty = 3)
abline(h = -2/sqrt(n), lty = 3)
it <- irregular.spec$rpspec * 10
ymin <- (min(irregular.spec$rpspec) - 1) * 10
ymax <- (max(irregular.spec$rpspec) + 1) * 10
plot(it, type = 'l',
main = "High Order AR-Spectrum as an approximation
to periodgram ( Order is fixed at 30 )",
ylab = "Rational Spectrum", xlab = "", ylim = c(ymin, ymax))
par(mfrow = c(1, 1))
# SPECTRUM OF DIFFERENCED ADJUSTED SERIES
lag <- length(adjusted.spec$acov) - 1
n <- adjusted.spec$n
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(3, 1))
plot((0:lag), adjusted.spec$acor, type = "h",
main="Autocorrelation & Parcor of Differenced Adjusted Series",
ylab = "Autocorrelation", xlab = "lag", ylim = c(-1, 1))
plot(adjusted.spec$parcor, type = "h", ylab = "Parcor", xlab = "lag",
ylim = c(-1, 1))
abline(h = 0, lty = 1)
abline(h = 2/sqrt(n), lty = 3)
abline(h = -2/sqrt(n), lty = 3)
it <- adjusted.spec$rpspec * 10
ymin <- (min(adjusted.spec$rpspec) - 1) * 10
ymax <- (max(adjusted.spec$rpspec) + 1) * 10
plot(it, type = 'l',
main = "High Order AR-Spectrum as an approximation
to periodgram( Order is fixed at 30 )",
xlab = "", ylab = "Rational Spectrum", ylim = c(ymin, ymax))
par(mfrow = c(1, 1))
# PARCOR OF TREND.ORDER TIME(S) DIFFERENCED TREND SERIES
lag <- length(differenced.trend$acov) - 1
n <- differenced.trend$n
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(2, 1))
plot((0:lag), differenced.trend$acor, type = "h",
main=paste(trend.order, "Time(s) Differenced Trend Series"),
ylab = "Autocorrelation", cex.main = 0.9, xlab = "lag",
ylim = c(-1, 1))
plot(differenced.trend$parcor, type = "h", ylab = "Parcor", xlab = "lag",
ylim = c(-1, 1))
abline(h = 0, lty = 1)
abline(h = 2/sqrt(n), lty = 3)
abline(h = -2/sqrt(n), lty = 3)
par(mfrow = c(1, 1))
# PARCOR OF SEASONAL.ORDER TIME(S) DIFFERENCED SEASONAL SERIES
lag <- length(differenced.season$acov) - 1
n <- differenced.season$n
if (separate.graphics == TRUE)
dev.new()
par(mfrow = c(2, 1))
plot((0:lag), differenced.season$acor, type = "h",
main=paste(seasonal.order, "time(s) Differenced Seasonal Series"),
cex.main = 0.9, ylab = "Autocorrelation", xlab = "lag",
ylim = c(-1, 1))
plot(differenced.season$parcor, type = "h", ylab = "Parcor",
xlab = "lag", ylim = c(-1, 1))
abline(h = 0, lty = 1)
abline(h = 2/sqrt(n), lty = 3)
abline(h = -2/sqrt(n), lty = 3)
}
par(oldpar)
}
decomp <- function(y, trend.order = 2, ar.order = 2, seasonal.order = 1,
period = 1, log = FALSE, trade = FALSE, diff = 1, miss = 0,
omax = 99999.9, plot = TRUE, ...)
{
if (trend.order < 1)
stop("trend.order is greater than or equal to 1")
m1 <- trend.order
m2 <- ar.order
ilog <- 0
if (log == TRUE)
ilog <- 1
itrade <- 0
if (trade == TRUE)
itrade <- 1
if (is.null(tsp(y)) == TRUE) {
year <- 1
month <- 1
} else if (is.null(tsp(y)) == FALSE) {
year <- start(y)[1]
month <- start(y)[2]
period <- tsp(y)[3]
}
if (is.null(tsp(y)[3]) == TRUE)
y <- ts(y, start = 1, frequency = period)
sorder <- seasonal.order
if (sorder != 0)
if (period < 2) {
warning("period must be greater than 1")
sorder <- 0
}
if (itrade != 0)
if(period != 4 && period !=12) {
warning("'trade = TRUE' is available only if period is 4 or 12")
itrade <- 0
}
n <- length(y)
ipar <- rep(0, 9)
ipar[1] <- trend.order
ipar[2] <- ar.order
ipar[3] <- period
ipar[4] <- sorder
ipar[5] <- ilog
ipar[6] <- itrade
ipar[7] <- diff
ipar[8] <- year
ipar[9] <- month
z <- .Fortran(C_decompf,
as.double(y),
as.integer(n),
as.integer(ipar),
trend = double(n),
seasonal = double(n),
ar = double(n),
trad = double(n),
noise = double(n),
para = double(26),
as.integer(miss),
as.double(omax),
ier = integer(1))
ier <- z$ier
if (ier != 0)
stop("log-transformation cannot be applied to zeros and nagative numbers")
trend <- z$trend
trend <- ts(trend, start = tsp(y)[1], frequency = tsp(y)[3])
season = z$seasonal
ar = z$ar
trading = z$trad
noise = z$noise
noise <- ts(noise, start = tsp(y)[1], frequency = tsp(y)[3])
aic = z$para[1]
lkhd = z$para[2]
sigma2 = z$para[3]
tau1 = z$para[4]
tau2 = z$para[5]
tau3 = z$para[6]
if (ar.order == 0 && sorder == 0) {
tau2 <- NULL
tau3 <- NULL
} else if (ar.order == 0 && sorder != 0) {
tau3 <- NULL
} else if (ar.order != 0 && sorder == 0) {
tau3 <- NULL
}
arcoef = z$para[7:(6+m2)]
tdf = z$para[(7+m2):(13+m2)]
if (sorder == 0)
season <- NULL
if (ar.order == 0)
ar <- NULL
if (itrade == 0)
trading <- NULL
conv.y <- y
if (log == TRUE)
conv.y <- log(conv.y)
if (miss > 0) {
for (i in 1:n)
if (y[i] > omax) {
conv.y[i] <- NA
noise[i] <- NA
}
} else if (miss < 0) {
for (i in 1:n)
if (y[i] < omax) {
conv.y[i] <- NA
noise[i] <- NA
}
}
out <- list(trend = trend, seasonal = season, ar = ar, trad = trading,
noise = noise, aic = aic, lkhd = lkhd, sigma2 = sigma2,
tau1 = tau1, tau2 = tau2, tau3 = tau3, arcoef = arcoef,
tdf = tdf, conv.y = conv.y)
class(out) <- "decomp"
if (plot == TRUE)
eval(parse(text=paste("plot.decomp(out, ...)")))
return(out)
}
#======================================================
plot.decomp <- function(x, ...)
#======================================================
{
y <- x$conv.y
ts.atr <- tsp(y)
trend <- x$trend
noise <- x$noise
n <- length(y)
period <- tsp(y)[3]
if (period == 4 || period == 12) {
xtime <- "year"
} else if (period == 24) {
xtime <- "day"
} else if (period == 5 || period == 7) {
xtime <- "week"
} else {
xtime <- "time"
}
oldpar <- par(no.readonly = TRUE)
nc <- 2
if (is.null(x$ar) == FALSE)
nc <- nc + 1
if (is.null(x$seasonal) == FALSE)
nc <- nc + 1
if (is.null(x$trad) == FALSE)
nc <- nc + 1
if (nc > 3)
par(mfrow = c((nc+1)/2, 2))
if (nc <= 3)
par(mfrow = c(nc, 1))
nmiss <- length(y) - length(y[!is.na(y)])
ymax <- max(y, na.rm = TRUE)
ymin <- min(y, na.rm = TRUE)
noise.nna <- noise[!is.na(noise)]
ymax <- max(trend, ymax)
ymin <- min(trend, ymin)
yd <- (ymax - ymin) / 10
ymax <- ymax + yd
ymin <- ymin - yd
### Original and Trend
if (nmiss == 0) {
plot(trend, type = 'l', col = 2, ylim = c(ymin, ymax),
main = "Original and Trend", xlab = "", ylab = "", ...)
par(new=TRUE)
plot(y, type = "l", ylim = c(ymin, ymax),
main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
} else {
yy1 <- rep(NA, n)
ymiss <- rep(NA, n)
ymiss <- ts(ymiss, start = ts.atr[1], frequency = ts.atr[3])
yy1 <- ts(yy1, start = ts.atr[1], frequency = ts.atr[3])
for (i in 1:n)
if (is.na(y[i]) == TRUE) {
ymiss[i] <- ymin
}
if ((is.na(y[1]) == FALSE) && (is.na(y[2]) = TRUE))
yy1[1] <- y[1]
if ((is.na(y[n-1]) == TRUE) && (is.na(y[n]) == FALSE))
yy1[n] <- y[n]
for(i in 2:(n-1)) {
if (is.na(y[i]) == FALSE)
if ((is.na(y[i-1]) == TRUE) && (is.na(y[i+1]) == TRUE))
yy1[i] <- y[i]
}
ymax <- ymax + 0.5 * yd
ymin <- ymin - 0.5 * yd
ymiss <- ymiss - 0.5 * yd
plot(trend, type = 'l', col = 2, ylim = c(ymin, ymax),
main = "Original and Trend", xlab = "", ylab = "", ...)
par(new=TRUE)
plot(y, type = "l", ylim = c(ymin, ymax),
main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
if (nmiss != 0) {
par(new = TRUE)
plot(yy1, type = "p", pch = 20, ylim = c(ymin, ymax), main = "",
xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex=0.6, ...)
par(new = TRUE)
plot(ymiss, type="p", pch = 20, col = 4, ylim = c(ymin, ymax),
main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n", cex=0.6,
...)
legend("topright", legend = "missing", col = 4, pch=20, bty = "n")
}
}
### Seasonal
ymax <- max(x$ar, x$seasonal, x$trad, noise.nna)
ymin <- min(x$ar, x$seasonal, x$trad, noise.nna)
my <- max(ymax, abs(ymin)) * 1.5
if (is.null(x$seasonal) == FALSE) {
season <- x$seasonal
season <- ts(season, start = ts.atr[1], frequency = ts.atr[3])
plot(season, type = "l", main = "Seasonal", xlab = "", ylab = "",
ylim = c(-my, my), ...)
}
### Noise
if (nmiss != 0) {
noise1 <- rep(NA, n)
noise1 <- ts(noise1, start = ts.atr[1], frequency = ts.atr[3])
if ((is.na(noise[1]) == FALSE) && (is.na(noise[2]) == TRUE))
noise1[1] <- noise[1]
if ((is.na(noise[n-1]) == TRUE) && (is.na(noise[n]) == FALSE))
noise1[n] <- noise[n]
for(i in 2:(n-1)) {
if (is.na(noise[i]) == FALSE)
if ((is.na(noise[i-1]) == TRUE) && (is.na(noise[i+1]) == TRUE))
noise1[i] <- noise[i]
}
plot(noise1, type = "p", pch = 20, main = "", xlab = "", ylab = "",
ylim = c(-my, my), cex=0.6, ...)
par(new = TRUE)
}
plot(noise, type = "l", main = "Noise", xlab = "", ylab = "",
ylim = c(-my, my), ...)
### AR
if (is.null(x$ar) == FALSE) {
ar <- x$ar
ar <- ts(ar, start = ts.atr[1], frequency = ts.atr[3])
plot(ar, type = "l", main = "AR component", xlab = "", ylab = "",
ylim = c(-my, my), ...)
}
### Trading Day Effect
if (is.null(x$trad) == FALSE) {
trad <- x$trad
trad <- ts(trad, start = ts.atr[1], frequency = ts.atr[3])
plot(trad, type = "l", main = "Trading Day Effect", xlab = "",
ylab = "", ylim = c(-my, my), ...)
}
par(mfrow=oldpar$mfrow, new = oldpar$new)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.