I am attempting to write an R function to aid in time series diagnostics.
The tsdiag() works well, but I would prefer to get a plot with ACF, PACF,
and Ljung-Box p-values of residuals. In my attempt to get to that point, I
am writing a function to calculate Ljung-Box p-values for residuals at
various lag distances.
ljung<-function(fit)
for(i in c(1:24,36,48))
{box<-(Box.test(fit$residuals, lag=i, type="Ljung-Box"))
print(c(i, box$p.value))}
This is one of my first R function writing attempts, so my apologies if
there is an obvious mistake. The above function produces the desired effect
in printing the lags and p-values to be plotted [where fit is the result of
arima()]; however I cannot seem to get the output stored in a data.frame for
subsequent plotting. I have tried storing the output using various methods
including data.frame, write.table, cat, capture.output, all with no success.
e.g:
ljung.out<-capture.output(print(c(i, box$p.value)))
#I saw a suggestion similar to this one in a previous post to the list...
Any hints on how I can get the output stored so that I may plot it later? I
am using v 1.9.1, RGui for Windows.
Thanks,
Andrew Kniss
Assistant Research Scientist
University of Wyoming
Dept. 3354
1000 E. University Ave.
Laramie, WY 82071
(307) 766-3949
akniss at uwyo.edu
Below is an arima fit that I have used in testing the function.
> dump("coal.fit", file=stdout())
coal.fit <-
structure(list(coef = structure(c(0.69539198614687, 484.903589344614
), .Names = c("ar1", "intercept")), sigma2 =
3109.45476265185,
var.coef = structure(c(0.0104024111201598, 0.302125085158484,
0.302125085158484, 634.39981868771), .Dim = as.integer(c(2,
2)), .Dimnames = list(c("ar1", "intercept"),
c("ar1", "intercept"
))), mask = c(TRUE, TRUE), loglik = -266.892361376605, aic 539.784722753209,
arma = as.integer(c(1, 0, 0, 0, 1, 0, 0)), residuals
structure(c(60.4342567589239,
-127.383559378086, -14.9885854976147, 123.839062585504,
-56.6019914334983,
35.7247594443981, 63.6906479431108, -28.1651273226733,
-6.91856808459544,
8.90309567990135, -30.8784722646861, -91.1489687772519,
-103.345257968621,
-29.2770349660464, -20.9664426335713, -25.3512422872431,
32.6086618928476, -6.98260117899269, -108.850345082021,
4.60267757422564,
38.6146462114696, 42.7187751257762, 79.9491758184327, 36.880952815858,
62.0132089128299, -0.848550671576206, -15.6420872534077,
111.955160137055, 13.5021374808082, -126.940710948639, 63.7127908071542,
27.4722158876983, -52.0448398629454, -15.4535767911051,
-73.4996569296364,
46.7008221699102, 27.5464232088949, -2.40151233395178,
-80.5337684309237,
-20.8162335807335, -18.2070175530272, -33.9885854976147,
-5.94848967770536, 17.8390625855041, 0.109559098069901,
39.5464232088949,
30.2537838322858, 32.9551601370546, 13.4381043864110), .Tsp = c(1,
49, 1), class = "ts"), call = stats:::arima(x = x, order = order,
seasonal = seasonal, include.mean = include.mean), series =
"x",
code = as.integer(0), n.cond = 0, model = structure(list(
phi = 0.69539198614687, theta = numeric(0), Delta = numeric(0),
Z = 1, a = 60.0964106553856, P = structure(0, .Dim = as.integer(c(1,
1))), T = structure(0.69539198614687, .Dim = as.integer(c(1,
1))), V = structure(1, .Dim = as.integer(c(1, 1))), h = 0,
Pn = structure(1, .Dim = as.integer(c(1, 1)))), .Names =
c("phi",
"theta", "Delta", "Z", "a",
"P", "T", "V", "h", "Pn")), x
structure(as.integer(c(569,
416, 422, 565, 484, 520, 573, 518, 501, 505, 468, 382, 310,
334, 359, 372, 439, 446, 349, 395, 461, 511, 583, 590, 620,
578, 534, 631, 600, 438, 516, 534, 467, 457, 392, 467, 500,
493, 410, 412, 416, 403, 422, 459, 467, 512, 534, 552, 545
)), .Dim = as.integer(c(49, 1)), .Dimnames = list(NULL, "Coal"),
.Tsp c(1,
49, 1), class = "ts")), .Names = c("coef",
"sigma2", "var.coef",
"mask", "loglik", "aic", "arma",
"residuals", "call", "series",
"code", "n.cond", "model", "x"), class =
"Arima")
Andrew Kniss wrote:> I am attempting to write an R function to aid in time series diagnostics. > The tsdiag() works well, but I would prefer to get a plot with ACF, PACF, > and Ljung-Box p-values of residuals. In my attempt to get to that point, I > am writing a function to calculate Ljung-Box p-values for residuals at > various lag distances. > > ljung<-function(fit) > for(i in c(1:24,36,48)) > {box<-(Box.test(fit$residuals, lag=i, type="Ljung-Box")) > print(c(i, box$p.value))} >You need to return() rather than print() the object. Uwe Ligges> This is one of my first R function writing attempts, so my apologies if > there is an obvious mistake. The above function produces the desired effect > in printing the lags and p-values to be plotted [where fit is the result of > arima()]; however I cannot seem to get the output stored in a data.frame for > subsequent plotting. I have tried storing the output using various methods > including data.frame, write.table, cat, capture.output, all with no success. > > > e.g: > ljung.out<-capture.output(print(c(i, box$p.value))) > #I saw a suggestion similar to this one in a previous post to the list... > > Any hints on how I can get the output stored so that I may plot it later? I > am using v 1.9.1, RGui for Windows. > Thanks, > > Andrew Kniss > Assistant Research Scientist > University of Wyoming > Dept. 3354 > 1000 E. University Ave. > Laramie, WY 82071 > (307) 766-3949 > akniss at uwyo.edu > > > Below is an arima fit that I have used in testing the function. > > >>dump("coal.fit", file=stdout()) > > coal.fit <- > structure(list(coef = structure(c(0.69539198614687, 484.903589344614 > ), .Names = c("ar1", "intercept")), sigma2 = 3109.45476265185, > var.coef = structure(c(0.0104024111201598, 0.302125085158484, > 0.302125085158484, 634.39981868771), .Dim = as.integer(c(2, > 2)), .Dimnames = list(c("ar1", "intercept"), c("ar1", "intercept" > ))), mask = c(TRUE, TRUE), loglik = -266.892361376605, aic > 539.784722753209, > arma = as.integer(c(1, 0, 0, 0, 1, 0, 0)), residuals > structure(c(60.4342567589239, > -127.383559378086, -14.9885854976147, 123.839062585504, > -56.6019914334983, > 35.7247594443981, 63.6906479431108, -28.1651273226733, > -6.91856808459544, > 8.90309567990135, -30.8784722646861, -91.1489687772519, > -103.345257968621, > -29.2770349660464, -20.9664426335713, -25.3512422872431, > 32.6086618928476, -6.98260117899269, -108.850345082021, > 4.60267757422564, > 38.6146462114696, 42.7187751257762, 79.9491758184327, 36.880952815858, > 62.0132089128299, -0.848550671576206, -15.6420872534077, > 111.955160137055, 13.5021374808082, -126.940710948639, 63.7127908071542, > > 27.4722158876983, -52.0448398629454, -15.4535767911051, > -73.4996569296364, > 46.7008221699102, 27.5464232088949, -2.40151233395178, > -80.5337684309237, > -20.8162335807335, -18.2070175530272, -33.9885854976147, > -5.94848967770536, 17.8390625855041, 0.109559098069901, > 39.5464232088949, > 30.2537838322858, 32.9551601370546, 13.4381043864110), .Tsp = c(1, > 49, 1), class = "ts"), call = stats:::arima(x = x, order = order, > seasonal = seasonal, include.mean = include.mean), series = "x", > code = as.integer(0), n.cond = 0, model = structure(list( > phi = 0.69539198614687, theta = numeric(0), Delta = numeric(0), > Z = 1, a = 60.0964106553856, P = structure(0, .Dim = as.integer(c(1, > > 1))), T = structure(0.69539198614687, .Dim = as.integer(c(1, > 1))), V = structure(1, .Dim = as.integer(c(1, 1))), h = 0, > Pn = structure(1, .Dim = as.integer(c(1, 1)))), .Names = c("phi", > "theta", "Delta", "Z", "a", "P", "T", "V", "h", "Pn")), x > structure(as.integer(c(569, > 416, 422, 565, 484, 520, 573, 518, 501, 505, 468, 382, 310, > 334, 359, 372, 439, 446, 349, 395, 461, 511, 583, 590, 620, > 578, 534, 631, 600, 438, 516, 534, 467, 457, 392, 467, 500, > 493, 410, 412, 416, 403, 422, 459, 467, 512, 534, 552, 545 > )), .Dim = as.integer(c(49, 1)), .Dimnames = list(NULL, "Coal"), .Tsp > c(1, > 49, 1), class = "ts")), .Names = c("coef", "sigma2", "var.coef", > "mask", "loglik", "aic", "arma", "residuals", "call", "series", > "code", "n.cond", "model", "x"), class = "Arima") > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> From: Uwe Ligges > > Andrew Kniss wrote: > > > I am attempting to write an R function to aid in time > series diagnostics. > > The tsdiag() works well, but I would prefer to get a plot > with ACF, PACF, > > and Ljung-Box p-values of residuals. In my attempt to get > to that point, I > > am writing a function to calculate Ljung-Box p-values for > residuals at > > various lag distances. > > > > ljung<-function(fit) > > for(i in c(1:24,36,48)) > > {box<-(Box.test(fit$residuals, lag=i, type="Ljung-Box")) > > print(c(i, box$p.value))} > > > > You need to return() rather than print() the object.Yes, but print() is supposed to return (invisibly) its argument(s)... Andy> Uwe Ligges > > > > This is one of my first R function writing attempts, so my > apologies if > > there is an obvious mistake. The above function produces > the desired effect > > in printing the lags and p-values to be plotted [where fit > is the result of > > arima()]; however I cannot seem to get the output stored in > a data.frame for > > subsequent plotting. I have tried storing the output using > various methods > > including data.frame, write.table, cat, capture.output, all > with no success. > > > > > > e.g: > > ljung.out<-capture.output(print(c(i, box$p.value))) > > #I saw a suggestion similar to this one in a previous post > to the list... > > > > Any hints on how I can get the output stored so that I may > plot it later? I > > am using v 1.9.1, RGui for Windows. > > Thanks, > > > > Andrew Kniss > > Assistant Research Scientist > > University of Wyoming > > Dept. 3354 > > 1000 E. University Ave. > > Laramie, WY 82071 > > (307) 766-3949 > > akniss at uwyo.edu > > > > > > Below is an arima fit that I have used in testing the function. > > > > > >>dump("coal.fit", file=stdout()) > > > > coal.fit <- > > structure(list(coef = structure(c(0.69539198614687, 484.903589344614 > > ), .Names = c("ar1", "intercept")), sigma2 = 3109.45476265185, > > var.coef = structure(c(0.0104024111201598, 0.302125085158484, > > 0.302125085158484, 634.39981868771), .Dim = as.integer(c(2, > > 2)), .Dimnames = list(c("ar1", "intercept"), c("ar1", > "intercept" > > ))), mask = c(TRUE, TRUE), loglik = -266.892361376605, aic > > 539.784722753209, > > arma = as.integer(c(1, 0, 0, 0, 1, 0, 0)), residuals > > structure(c(60.4342567589239, > > -127.383559378086, -14.9885854976147, 123.839062585504, > > -56.6019914334983, > > 35.7247594443981, 63.6906479431108, -28.1651273226733, > > -6.91856808459544, > > 8.90309567990135, -30.8784722646861, -91.1489687772519, > > -103.345257968621, > > -29.2770349660464, -20.9664426335713, -25.3512422872431, > > 32.6086618928476, -6.98260117899269, -108.850345082021, > > 4.60267757422564, > > 38.6146462114696, 42.7187751257762, 79.9491758184327, > 36.880952815858, > > 62.0132089128299, -0.848550671576206, -15.6420872534077, > > 111.955160137055, 13.5021374808082, -126.940710948639, > 63.7127908071542, > > > > 27.4722158876983, -52.0448398629454, -15.4535767911051, > > -73.4996569296364, > > 46.7008221699102, 27.5464232088949, -2.40151233395178, > > -80.5337684309237, > > -20.8162335807335, -18.2070175530272, -33.9885854976147, > > -5.94848967770536, 17.8390625855041, 0.109559098069901, > > 39.5464232088949, > > 30.2537838322858, 32.9551601370546, 13.4381043864110), > .Tsp = c(1, > > 49, 1), class = "ts"), call = stats:::arima(x = x, > order = order, > > seasonal = seasonal, include.mean = include.mean), > series = "x", > > code = as.integer(0), n.cond = 0, model = structure(list( > > phi = 0.69539198614687, theta = numeric(0), Delta = > numeric(0), > > Z = 1, a = 60.0964106553856, P = structure(0, .Dim > = as.integer(c(1, > > > > 1))), T = structure(0.69539198614687, .Dim = > as.integer(c(1, > > 1))), V = structure(1, .Dim = as.integer(c(1, 1))), h = 0, > > Pn = structure(1, .Dim = as.integer(c(1, 1)))), > .Names = c("phi", > > "theta", "Delta", "Z", "a", "P", "T", "V", "h", "Pn")), x > > structure(as.integer(c(569, > > 416, 422, 565, 484, 520, 573, 518, 501, 505, 468, 382, 310, > > 334, 359, 372, 439, 446, 349, 395, 461, 511, 583, 590, 620, > > 578, 534, 631, 600, 438, 516, 534, 467, 457, 392, 467, 500, > > 493, 410, 412, 416, 403, 422, 459, 467, 512, 534, 552, 545 > > )), .Dim = as.integer(c(49, 1)), .Dimnames = list(NULL, > "Coal"), .Tsp > > c(1, > > 49, 1), class = "ts")), .Names = c("coef", "sigma2", > "var.coef", > > "mask", "loglik", "aic", "arma", "residuals", "call", "series", > > "code", "n.cond", "model", "x"), class = "Arima") > > > > ______________________________________________ > > R-help at stat.math.ethz.ch mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > >