R version 2.14.1 (2011-12-22) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. [Previously saved workspace restored] > #Run from command line using R CMD BATCH ar1_example.r > #Run from inside R using source("ar1_example.r",echo=T) > > num <- 100 > desired_width <- 800 > desired_height <- 800 > desired_res <- 120 > > legends <- c("phi=0.0", "phi=0.2", "phi=0.9","Trend + 'phi=0.2'") > pt_two_col <- "forestgreen" > colors <- c("black", pt_two_col, "red", pt_two_col) > width <- 3 > widths <- c(width,width,width,width) > ltys <- c(2,1,6,3) > actualvalues <- c(0.0,0.2,0.9,0.2) > white <- arima.sim(n = num, list(ar = 0.0), innov=rnorm(num)) Warning message: In min(Mod(polyroot(c(1, -model$ar)))) : no non-missing arguments to min; returning Inf > pt_two <- arima.sim(n = num, list(ar = 0.2), innov=rnorm(num)) > pt_nine <- arima.sim(n = num, list(ar = 0.9), innov=rnorm(num)) > signal_plus_pt_two <- pt_two + 0.1*(1:num) > white <- white - mean(white) > pt_two <- pt_two - mean(pt_two) > pt_nine <- pt_nine - mean(pt_nine) > signal_plus_pt_two <- signal_plus_pt_two - mean(signal_plus_pt_two) > png(file="fig1.png",width=desired_width,height=desired_height,res=desired_res) > ts.plot(white, pt_two, pt_nine, signal_plus_pt_two, gpars=list(xlab="time", ylab="", main="AR(1) noise time series with different parameters 'phi'", col=colors, lty=ltys, lwd = widths)) > legend("topleft", legend = legends, col=colors, lty=ltys, lwd=widths, cex=1) > dev.off function (which = dev.cur()) { if (which == 1) stop("cannot shut down device 1 (the null device)") .Internal(dev.off(as.integer(which))) dev.cur() } > > print(.Random.seed) [1] 403 208 -1486295240 -576719525 586825717 -1994543133 [7] -699128604 -872001622 597633127 -1089920975 -1766573397 -725153155 [13] -708710169 -1582122374 381692341 -292651968 -651968843 -265504729 [19] -1060241936 -1535742199 906109086 -1971457012 -514321940 -1340053167 [25] -1976868459 -2069502369 130803388 97337877 -795450146 788684202 [31] 1760560427 529389170 654228630 164741407 1930167317 1302540105 [37] 1444937816 -393620445 2036606688 -1920828539 890501186 694739334 [43] 1502340510 1601008347 -726268809 1433958273 -708507144 134597919 [49] -880231236 -1829109376 -1404964349 1472981884 1175889452 1199315972 [55] -1586641248 -2034221892 -2050859357 1734563007 519511373 1923179800 [61] -2047722002 778450710 1470292036 217742795 1267417512 1582803173 [67] 215205574 -1314376201 1363266979 1134227416 1803890479 -59478145 [73] -959731171 1341131187 1167331516 -1896009288 1613661355 1376171293 [79] 732290059 1910529936 34122202 85750338 -1306436684 -857233418 [85] 1219426937 208515114 -1309463535 1250152492 -597591260 -685786490 [91] 1844171005 -1912775419 -417597332 839484945 224833765 1539871632 [97] -59200612 -305398845 41037313 1855803630 1077868340 1914825377 [103] -390551904 576570413 -639175689 -1337024344 -696071222 1156759572 [109] 1867220107 728720551 1072066685 1644275888 -1645312022 -1899839982 [115] -828339649 1099218294 -878793614 -1043306889 -1376182963 -1656195960 [121] -2075757812 415226219 -1262247457 1685911567 -735523925 -770105868 [127] -96328230 988300656 -1471315702 195301749 840910629 154769509 [133] -1428629511 1517610858 -1313108691 159974569 627268576 -977533071 [139] -863999595 -1457583863 -1574669518 761501654 -1718726222 -1783470337 [145] -619622743 582259413 -751012514 327987033 1424518776 -1723712296 [151] 185676615 -44952673 430797259 253695587 1258867368 -1874038822 [157] 469495887 1017472464 1803765963 -941430507 2019399016 1463547656 [163] 1483762191 389111872 61685697 -465938773 -882923949 -130317099 [169] 1567477587 -1182848588 711567734 445724564 -508746334 1048501038 [175] -1256329930 959543376 140345719 -912858396 788283856 35951828 [181] 1387905478 621818998 1595912595 -1981735653 -1553389520 1463360668 [187] 782033816 2083856794 -742778271 1363405567 -1629922334 -1891114920 [193] -1058713347 1468015398 -1002858558 1032501924 1559878079 646823804 [199] 896661633 -139704376 614142276 -7905462 856197501 -948294376 [205] -1631459967 -273890201 1900575419 -19791396 390346769 1947642540 [211] -1621232875 1379450885 -1010624069 -502724659 2043250729 1556759529 [217] 490091692 -807857120 -2092400996 -873032499 -450879542 684081637 [223] -1037230907 1262907353 2046851638 -1943919412 -832981512 1919730444 [229] -1181995659 1729323849 -1386398185 1562300524 -868495568 925145895 [235] -1547461075 754067700 1040977453 -1013511951 -1769784445 1855845040 [241] -761622780 1123431920 971192514 1107553998 -1706452099 177824715 [247] 1430455293 349598091 -1827385706 816500414 765310528 -171844524 [253] 1040208548 434694517 1595325151 -1243435940 421981146 -172914788 [259] -1325185940 2128101715 1952928335 -966587028 -9910401 -20440776 [265] -1317710292 -1018381603 2039935430 -269932997 -254608217 -1338450162 [271] -664880735 -220523228 927341116 2106172989 655138566 -589831395 [277] 551169340 -808108895 2113002472 -396236331 394469109 1672515852 [283] -1586512200 243412780 1705809967 483633709 1203193170 1432196136 [289] 2047499560 -1563895777 1522344842 1154144262 -1421203059 -1737633843 [295] 1535227010 -1472880831 808233981 -1303779064 591043 1208619727 [301] -380369358 -1380368620 1416776038 1176059716 864962386 503019802 [307] -1350333312 -684531931 1055417590 -1943454667 1890134699 1700585433 [313] 265707779 -1619167589 1949626311 -1086685579 811428918 1476164251 [319] -1950355322 591766935 178806726 1809978698 1509688141 2136513455 [325] -748404326 1474035097 -970733115 1525582760 -1978957920 825055067 [331] 153401975 757716383 -940304298 -1101698694 -1351251153 10423686 [337] 1510712592 -1376240133 1200579608 -464160792 31880545 1928150728 [343] 788559691 -1509427206 618973720 -1361677025 -1867489274 39104121 [349] -596758942 1330718121 -990810514 -1616429666 -847157795 -1244913679 [355] 1922146835 -1439476545 1767826671 1448625005 1797069134 -1130977850 [361] -2050165081 -1569085272 -570611074 -1716837046 -2118631735 -737158219 [367] -578744610 -469590147 -1088640839 1883879047 -1829841769 -357127222 [373] 628384540 -1181929936 1360186325 1840832877 -1184556007 417042250 [379] -867738642 -423430784 468360083 -877022816 2124390792 -218589825 [385] -1360259137 1298209587 2054757557 1594220533 1736722073 -2057032089 [391] 1763638937 -1722085682 700541178 141930628 -1792050636 -343366956 [397] -182564464 -912163920 704871418 137718421 888129513 -688526420 [403] -787417568 -803091476 -1094998152 -1134453177 90764850 1284413666 [409] 1225346528 -1926093314 1315333749 -910228880 -1376681403 -1027836235 [415] 937186007 -1463182613 1943042880 1091441109 -408766266 1095370245 [421] 1622001203 -1920294803 746849459 1421876696 63639852 -109136225 [427] -77523893 1497805908 1088617470 -696922867 552122424 964649866 [433] -1689255801 -1979765355 -2067614983 -648292530 -156908935 186675731 [439] -2047249362 -1222427412 717319385 -278592334 -434121187 -934187001 [445] -1734893625 -946659615 -1765185150 -577828832 -1366434599 864246406 [451] 1786345703 1282065372 -1502405071 32966866 8327348 225233678 [457] 1695293629 118438805 -557558371 -1658422976 15694480 -1389754332 [463] -843515024 -1679962788 1302613531 -168456913 -988187061 1337818020 [469] 2040019126 -326882474 -1162386354 -1576768674 -1121852119 -441973051 [475] 1634664385 983359761 405097117 -286122899 131515818 1757092988 [481] -1599519466 -532396572 -1837702162 1247191513 -1756960198 1425988491 [487] -1033395177 781623300 2127364887 1167255330 -1626259899 1948878197 [493] 978016279 -74726965 -2085275517 666550449 1219786296 -1197479121 [499] -1902807933 1765942445 1180911171 -338979497 -1314557351 -1639736194 [505] -1150793527 -426477575 1575216958 -617595158 61255303 1839301945 [511] 1726054704 1468185480 -423052619 1613644939 -282998276 1911870047 [517] -1705842719 1946828139 -1599058489 88267733 234269841 1264843688 [523] -291220214 1839558239 283791452 230859162 1883848528 90413516 [529] 1266155734 2072971476 2124307697 1474225968 -598957390 -1012463268 [535] -1968137285 1169992235 -1596835365 -1239618899 1978324809 -693062490 [541] -1639393424 231913531 -617169071 -716598064 1918233110 -1435084971 [547] -1787068827 -1723880944 1569987643 -505313652 507739183 187311426 [553] 1016666753 689756089 -993108049 -832950202 -1573743722 1760883187 [559] -1800487977 2056205568 574497289 2000738271 458883846 21411408 [565] 1912206441 -1708703237 -1800443530 -1402167844 -222218455 -1780156312 [571] 950413766 -1599081957 1306519266 -904607010 815348545 2126166973 [577] -1936923655 -931396869 -1601149931 -1868820323 580727686 -1616216927 [583] 1210718496 160574491 1862149934 -1728414535 -44476584 1164321624 [589] 2090171424 -1668620131 914531256 -1978272510 1388267452 -892128775 [595] 1634998612 -759564644 -1906490239 371422694 -652180123 -1499564930 [601] -869957364 255997405 -1525799062 -756369516 540282492 498145101 [607] 845041190 -987407192 1945458898 84293929 532601702 -224044446 [613] -22727092 -1775749221 -502906449 552539214 1843608398 -627680460 [619] -1651793254 -1657632360 -1608882038 -1907956082 -1391233504 1444431158 [625] -1294523550 968297574 > > white_noise <- arima(white,order=c(1,0,0)) > pt_two_noise <- arima(pt_two,order=c(1,0,0)) > pt_nine_noise <- arima(pt_nine,order=c(1,0,0)) > signal_plus_noise <- arima(signal_plus_pt_two,order=c(1,0,0)) > print(white_noise) Call: arima(x = white, order = c(1, 0, 0)) Coefficients: ar1 intercept 0.0890 0.0009 s.e. 0.0993 0.1014 sigma^2 estimated as 0.8553: log likelihood = -134.08, aic = 274.16 > print(pt_two_noise) Call: arima(x = pt_two, order = c(1, 0, 0)) Coefficients: ar1 intercept 0.1830 -0.0006 s.e. 0.0978 0.1353 sigma^2 estimated as 1.228: log likelihood = -152.18, aic = 310.36 > print(pt_nine_noise) Call: arima(x = pt_nine, order = c(1, 0, 0)) Coefficients: ar1 intercept 0.891 0.2011 s.e. 0.049 0.9147 sigma^2 estimated as 1.145: log likelihood = -149.46, aic = 304.92 > print(signal_plus_noise) Call: arima(x = signal_plus_pt_two, order = c(1, 0, 0)) Coefficients: ar1 intercept 0.9035 -0.0205 s.e. 0.0434 1.3483 sigma^2 estimated as 2.01: log likelihood = -177.64, aic = 361.28 > > png(file="fig2.png",width=desired_width,height=desired_height,res=desired_res) > plot(acf(white,plot=FALSE), type="l", main="Autocorrelation function (ACF) of AR(1) time series", max.mfrow=1,col=colors[1], xlab="lag", ylab="", lty=ltys[1], lwd = widths[1]) > lines(acf(pt_two,plot=FALSE)$acf[-1,1,1], lty=ltys[2], col=colors[2], lwd=widths[2]) > lines(acf(pt_nine,plot=FALSE)$acf[-1,1,1], lty=ltys[3], col=colors[3], lwd=widths[3]) > lines(acf(signal_plus_pt_two,plot=FALSE)$acf[-1,1,1], lty=ltys[4], col=colors[4], lwd=widths[4]) > legend("topright", legend = legends, col=colors, lty=ltys, lwd=widths, cex=1) > dev.off function (which = dev.cur()) { if (which == 1) stop("cannot shut down device 1 (the null device)") .Internal(dev.off(as.integer(which))) dev.cur() } > > means <- c(white_noise$coef[1], pt_two_noise$coef[1], pt_nine_noise$coef[1], signal_plus_noise$coef[1]) > names <- legends > standardErrors <- c(sqrt(diag(vcov(white_noise))) [1], sqrt(diag(vcov(pt_two_noise))) [1], sqrt(diag(vcov(pt_nine_noise))) [1], sqrt(diag(vcov(signal_plus_noise))) [1]) > > #Plot actual and estimated AR(1) parameters. > offset <- 0.2 > width <- 1 > separation <- 1.2 > actualhalfwidth <- 0.02 > lowest <- min(c(0,min(means - 2*standardErrors))) - offset/2 > highest <- max(means + 2*standardErrors) > plotx <- offset + c(0,width,width,0,NA,separation,separation+width,separation+width,separation,NA,2*separation,2*separation+width,2*separation+width,2*separation,NA,3*separation,3*separation+width,3*separation+width,3*separation) > ploty1 <- c(means[1]+2*standardErrors[1],means[1]+2*standardErrors[1],means[1]-2*standardErrors[1],means[1]-2*standardErrors[1],NA,means[2]+2*standardErrors[2],means[2]+2*standardErrors[2],means[2]-2*standardErrors[2],means[2]-2*standardErrors[2],NA,means[3]+2*standardErrors[3],means[3]+2*standardErrors[3],means[3]-2*standardErrors[3],means[3]-2*standardErrors[3],NA,means[4]+2*standardErrors[4],means[4]+2*standardErrors[4],means[4]-2*standardErrors[4],means[4]-2*standardErrors[4]) > ploty2 <- c(actualvalues[1]+actualhalfwidth,actualvalues[1]+actualhalfwidth,actualvalues[1]-actualhalfwidth,actualvalues[1]-actualhalfwidth,NA,actualvalues[2]+actualhalfwidth,actualvalues[2]+actualhalfwidth,actualvalues[2]-actualhalfwidth,actualvalues[2]-actualhalfwidth,NA,actualvalues[3]+actualhalfwidth,actualvalues[3]+actualhalfwidth,actualvalues[3]-actualhalfwidth,actualvalues[3]-actualhalfwidth,NA,actualvalues[4]+actualhalfwidth,actualvalues[4]+actualhalfwidth,actualvalues[4]-actualhalfwidth,actualvalues[4]-actualhalfwidth) > png(file="fig3.png",width=desired_width,height=desired_height,res=desired_res) > barCenters <- barplot(means, names.arg=names, cex.names=1.2, yaxt="n", border=NA, col=NA, las=1, ylim=c(lowest,highest)) > axis(side = 2, at = c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0), las=1) > title(main = "Actual vs estimated AR(1) parameters 'phi'") > polygon(plotx, ploty1, density = NA, border = "grey", col = "grey") > polygon(plotx, ploty2, density = NA, border = colors, col = colors) > dev.off function (which = dev.cur()) { if (which == 1) stop("cannot shut down device 1 (the null device)") .Internal(dev.off(as.integer(which))) dev.cur() } > > proc.time() user system elapsed 0.488 0.040 0.521