Chatfield’s Plots in S-Plus

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..
1.2
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)])
About these ads
This entry was posted in R and tagged , , . Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s