Chatfield’s Plots in S-Plus

July 9, 2010
By

(This article was first published on Guy Abel » R, and kindly contributed to R-bloggers)

I have recently finished reading the sixth edition of The Analysis of Time Series: An Introduction by Chatfield in our Statistics reading group. Whilst enjoying most of the book I got a little confused when looking at Appendix D: Some MINITAB and S-PLUS Commands. In the S-Plus section the author gives the code below to replicate his Figure 1.2.

postscript(file="1.2.ps")
recife<-scan("Recife")
RecfS<-cts(recife, start=dates("01/01/1953"),units="months",frequency="12")
fst_as.numeric(date("01/01/1953"),format="dd/mm/yyyy")
mid1_as.numeric(dates("01/01/1954"),format="dd/mm/yyyy")
mid2_as.numeric(dates("01/01/1955"),format="dd/mm/yyyy")
mid3_as.numeric(dates("01/01/1956"),format="dd/mm/yyyy")
mid4_as.numeric(dates("01/01/1957"),format="dd/mm/yyyy")
mid5_as.numeric(dates("01/01/1958"),format="dd/mm/yyyy")
mid6_as.numeric(dates("01/01/1959"),format="dd/mm/yyyy")
mid7_as.numeric(dates("01/01/1960"),format="dd/mm/yyyy")
mid8_as.numeric(dates("01/01/1961"),format="dd/mm/yyyy")
lst_as.numeric(dates("01/01/1962"),format="dd/mm/yyyy")
ts.plot(RecfS,xaxt="n",ylab="Temperature (deg C)", xlab="Year", type="l")
axis(side=1, at = c(fst,mid1,mid2,mid3,mid4,mid5,mid6,mid7,mid8,lst),
labels=c("Jan 53", "Jan 54", "Jan 55", "Jan 56", "Jan 57", "Jan 58",
"Jan 59", "Jan 60", "Jan 61", "Jan 62"),
ticks=T)
dev.off()

I was not too sure what was going on with the code, which gave a tonne of error messages, some of which might well have been typo’s by the publisher (_ instead of <-)? In addition, the author stresses the effort required to construct nice plots in S-Plus. This got me thinking that 1) his code was excessive (not just because it does not work) and 2) he could have saved a lot of his effort by using R. The R code below proves my point

> png("...:/1.2.png", width=650, height=460, units="px")
> recife <- scan("http://people.bath.ac.uk/mascc/Recife.TS", sep="", skip=3)
> n<-length(recife)
> times<-seq(as.Date("1953/1/1"), by="month", length.out=n)
> plot(recife, xaxt="n", type= "l",ylab="Temperature (deg C)", xlab="Year")
> axis(1, at = seq(1,n,12), labels = format(times, "%b %Y")[seq(1,n,12)])
> dev.off()

Much Tidier! Here is the plot..Note, if you only wanted labels every other January, as in p.2 of the book (but not in the S-Plus code), then you can use..

> axis(1, at = seq(1,n,24), labels = format(times, "%b %y")[seq(1,n,24)])

To leave a comment for the author, please follow the link and comment on his blog: Guy Abel » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.