{"url": "https://hombao.ics.uci.edu/R-Computing-TIMEDOMAIN.txt", "content": "\n## Install astsa (one time only)\n\n## Load astsa (every time you need to use this package)\n\n## To access all datasets in astsa\nastsadata()\n\n## These commands are taken from Shumway and Stoffer (2010)\n\n############### Global temperature data ####################\n\n## plot\npar(mfrow=c(1,1))\nplot(gtemp, type=\"o\", ylab=\"Global Temperature Deviations\")\n\n## fit a linear model\nfit <- lm(gtemp ~ time(gtemp)) ## regress gtemp on time\n## fit is an object that stores output from lm\nnames(fit)\nsummary(fit)\npar(mfrow=c(1,1))\nplot(gtemp, type=\"o\", ylab=\"Global Temperature Deviation\")\nabline(fit, col=2) ## add estimated regression line to the plot\n\n## examine the residuals\nresid = fit$resid;\npred = fit$fitted;\npar(mfrow=c(2,1)); \nplot(resid, type=\"o\", ylab=\"Residuals\", xlab=\"time\"); \nN = length(resid);\nzeroline = rep(0, N);\nlines(zeroline, col=2);\nplot(y=resid, x=pred, type=\"o\", ylab=\"Residuals\", xlab=\"predicted\"); \nlines(x=pred, y=zeroline, col=2);\n\npar(mfrow=c(2,1)); \nacf(resid);\npacf(resid); \n\n\n############### LA County Cardiac Mortality Data ##############\n\n## Plot the LA County mortality dataset\n## Series: cmort, temperature, particulate\npar(mfrow=c(3,1))\nplot(cmort, main=\"Cardiovascular Mortality\", xlab=\"\", ylab=\"\"); \nplot(tempr, main=\"Temperature\", xlab=\"\", ylab=\"\")\nplot(part, main=\"Particulates\", xlab=\"\", ylab=\"\")\n\ndev.new()\npairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part))\n\n##### Fit a model for cardiac mortality\ntemp= tempr-mean(tempr) ## center temperature \ntemp2 = temp^2        ## square of temp\ntrend = time(cmort)  ## time\nfit = lm(cmort~trend + temp + temp2 + part, na.action=NULL)\nsummary(fit)    ## regression results\nsummary(aov(fit))  ## ANOVA table   (compare to next line) \nsummary(aov(lm(cmort~cbind(trend, temp, temp2, part)))) ##  Table 2.1 \n\n## Illustration on how to compute the information criteria\nnum = length(cmort) ##  sample size\nAIC(fit)/num- log(2*pi) ## AIC\nAIC(fit, k=log(num))/num- log(2*pi) ##  BIC\n(AICc = log(sum(resid(fit)^2)/num)+ (num+5)/(num-5-2)) ## AICc\n\n## Example on smoothing/filtering\nma5 = filter(cmort, sides=2, rep(1,5)/5)\nma53 = filter(cmort, sides=2, rep(1,53)/53) \npar(mfrow=c(1,1));\nplot(cmort, type=\"p\", ylab=\"mortality\") \nlines(ma5,col=2)\nlines(ma53,col=3)\n\n##### Fitting some polynomial and sinusoidal trends\nwk = time(cmort) - mean(time(cmort)) ## wk is essentially t/52 centered at zero\nwk2 = wk^2 \nwk3 = wk^3\ncs = cos(2*pi*wk)\nsn = sin(2*pi*wk)\nreg1 = lm(cmort ~ wk + wk2 + wk3, na.action=NULL)\nreg2 = lm(cmort ~ wk + wk2 + wk3 + cs + sn, na.action=NULL) \npar(mfrow=c(1,1));\nplot(cmort, type=\"p\", ylab=\"mortality\") \nlines(fitted(reg1), col=2)\nlines(fitted(reg2), col=3)\n\n## Compare the residual series for each model\nresid1 = reg1$resid\nresid2 = reg2$resid\npar(mfrow=c(2,1));\nplot(resid1)\nplot(resid2)\n\npar(mfrow=c(2,2)); \nacf(resid1)\nacf(resid2)\npacf(resid1)\npacf(resid2)\n\n\n###### Kernel smoothing (with Guassian kernel)\npar(mfrow=c(1,1));\nplot(cmort, type=\"p\", ylab=\"mortality\");\nsm1 = ksmooth(time(cmort), cmort, \"normal\", bandwidth=5/52);\nsm2 = ksmooth(time(cmort), cmort, \"normal\", bandwidth=2);\nlines(sm1, col=2); \nlines(sm2, col=3);\n\n\n#################### ARIMA ##########################\n\n##### Generating time series \n\n##### AR\nts1 = arima.sim(list(order=c(1,0,0), ar=0.9), n=100);\nts2 = arima.sim(list(order=c(1,0,0), ar=-0.9), n=100);\n\npar(mfrow=c(2,1))  ## in the expression below, -is a  space and == is equal\nplot(ts1, ylab=\"TS\", xlab=\"Time\"); title(paste(\"AR(1) phi=+.9\"));\nplot(ts2, ylab=\"TS\", xlab=\"Time\"); title(paste(\"AR(1) phi=-.9\"));\n\n##### MA\nts1 = arima.sim(list(order=c(0,0,1), ma=c(0.5)), n=100);\nts2 = arima.sim(list(order=c(0,0,5), ma=c(1, 0.3, 0.2, 0.1, 0.05)), n=100);\n\npar(mfrow=c(2,1))  ## in the expression below, -is a  space and == is equal\nplot(ts1, ylab=\"TS\", xlab=\"Time\"); title(paste(\"MA(1)\"));\nplot(ts2, ylab=\"TS\", xlab=\"Time\"); title(paste(\"MA(5)\"));\n\n\n##### Fitting ARIMA models to data\n\n## Data is generated by an AR(1) process\ndata1 = arima.sim(list(order=c(1,0,0), ar=0.9), n=100);\n\n## Fit ARMA(2,0,1) to data1\nfit1 = arima(data1, order=c(2,0,1))\n\n\n\n\n\n\n\n\n\n", "encoding": "ascii"}