R version 2.8.0 (2008-10-20) Copyright (C) 2008 The R Foundation for Statistical Computing ISBN 3-900051-07-0 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. > x <- c(186448,190530,194207,190855,200779,204428,207617,212071,214239,215883,223484,221529,225247,226699,231406,232324,237192,236727,240698,240688,245283,243556,247826,245798,250479,249216,251896,247616,249994,246552,248771,247551,249745,245742,249019,245841,248771,244723,246878,246014,248496,244351,248016,246509,249426,247840,251035,250161,254278,250801,253985,249174,251287,247947,249992,243805,255812,250417,253033,248705,253950,251484,251093,245996,252721,248019,250464,245571,252690,250183,253639,254436,265280,268705,270643,271480) > par1 = '4' > #'GNU S' R Code compiled by R2WASP v. 1.0.44 () > #Author: Prof. Dr. P. Wessa > #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ > #Source of accompanying publication: Office for Research, Development, and Education > #Technical description: Write here your technical program description (don't use hard returns!) > par1 <- as.numeric(par1) > nx <- length(x) > x <- ts(x,frequency=par1) > m <- StructTS(x,type='BSM') > m$coef level slope seas epsilon 107560.8 631968.3 148365.8 1311757.2 > m$fitted level slope sea 1 Q1 186448.0 0.00000 0.00000 1 Q2 189443.2 1029.33220 890.02079 1 Q3 192177.6 1813.07691 1363.11207 1 Q4 191178.4 367.99751 726.73214 2 Q1 198822.7 3969.89692 -55.64647 2 Q2 203897.8 4508.18807 175.67714 2 Q3 207341.3 3983.15050 632.88766 2 Q4 212275.9 4455.45676 -520.64435 3 Q1 215377.1 3803.47773 -719.38139 3 Q2 216768.8 2621.69111 -135.93830 3 Q3 221608.5 3705.80559 1168.66381 3 Q4 223102.8 2619.75286 -864.77165 4 Q1 225552.6 2537.20478 -252.54159 4 Q2 227521.9 2259.83846 -645.36932 4 Q3 229793.6 2265.60041 1608.68688 4 Q4 232753.4 2604.87643 -648.50572 5 Q1 236619.8 3218.62545 178.31118 5 Q2 238311.7 2473.94315 -1106.97633 5 Q3 239785.9 1986.87614 1225.24823 5 Q4 241610.6 1907.76975 -871.68601 6 Q1 244306.9 2291.68662 729.78506 6 Q2 245259.7 1639.14749 -1285.03504 6 Q3 246673.1 1529.17384 1223.48491 6 Q4 247302.4 1090.58817 -1222.71708 7 Q1 249058.4 1414.63283 1212.76222 7 Q2 250462.9 1409.73124 -1243.78079 7 Q3 251009.2 989.17844 1156.60863 7 Q4 250057.0 43.33840 -1834.01173 8 Q1 249202.1 -394.15284 1072.53394 8 Q2 248002.3 -786.64065 -1198.46316 8 Q3 247188.6 -799.79933 1590.82185 8 Q4 248161.9 63.93423 -1165.11644 9 Q1 248447.6 171.95929 1228.05452 9 Q2 247586.5 -331.29186 -1521.60847 9 Q3 247569.5 -178.19104 1351.28216 9 Q4 247164.3 -288.76545 -1252.36481 10 Q1 247168.4 -146.12126 1511.09633 10 Q2 246571.4 -365.72048 -1707.56106 10 Q3 245740.9 -592.14601 1282.39686 10 Q4 246510.9 71.38755 -922.62481 11 Q1 246804.0 179.33370 1622.79683 11 Q2 246366.5 -121.09267 -1822.77902 11 Q3 246717.4 108.79679 1151.13567 11 Q4 247238.5 309.65253 -858.39872 12 Q1 247649.6 359.05532 1744.68297 12 Q2 249127.6 904.04898 -1637.18863 12 Q3 250017.7 897.26148 1021.67244 12 Q4 251013.3 945.16409 -883.01899 13 Q1 252467.0 1192.86115 1652.11541 13 Q2 252896.4 820.98123 -1856.79867 13 Q3 253242.6 589.76429 890.68001 13 Q4 251477.0 -557.51373 -1567.01868 14 Q1 249945.1 -1032.14657 1646.42452 14 Q2 249354.4 -817.08587 -1545.40113 14 Q3 248590.0 -791.46452 1385.61043 14 Q4 246144.3 -1597.17670 -1822.45975 15 Q1 250668.9 1384.61779 3230.25899 15 Q2 252114.6 1414.38936 -1716.73148 15 Q3 252233.9 783.55594 1203.78821 15 Q4 252184.3 377.72961 -3218.93511 16 Q1 251359.0 -208.25108 2966.94635 16 Q2 252245.1 324.80784 -1103.06677 16 Q3 250767.6 -553.05201 888.52228 16 Q4 249405.4 -947.17194 -3156.60624 17 Q1 249369.0 -503.58645 3067.47217 17 Q2 248806.8 -532.12546 -769.47751 17 Q3 249012.1 -172.93206 1221.47644 17 Q4 248899.7 -143.42893 -3347.66322 18 Q1 249312.8 127.62633 3203.32114 18 Q2 250507.2 647.21653 -657.48094 18 Q3 251997.9 1058.06169 1377.58289 18 Q4 256164.3 2572.11654 -2699.60969 19 Q1 261077.4 3712.37150 3471.08344 19 Q2 267906.4 5230.40448 -175.21512 19 Q3 271186.9 4280.65434 65.35888 19 Q4 274956.4 4031.67375 -3316.67904 > m$resid Qtr1 Qtr2 Qtr3 Qtr4 1 0.000000000 0.782514217 1.197272205 -1.834519909 2 4.366719796 0.675983271 -0.672827275 0.591445456 3 -0.837915512 -1.459492842 1.371302334 -1.364392979 4 -0.104653900 -0.346615297 0.007267340 0.426438434 5 0.774292655 -0.934477533 -0.613374916 -0.099471635 6 0.483466684 -0.820104954 -0.138401058 -0.551613513 7 0.407789111 -0.006163784 -0.529117614 -1.189708005 8 -0.550412173 -0.493658962 -0.016553720 1.086476898 9 0.135894325 -0.633021915 0.192593574 -0.139092154 10 0.179438389 -0.276233734 -0.284827692 0.834667888 11 0.135788418 -0.377909849 0.289183195 0.252659775 12 0.062144819 0.685556227 -0.008538124 0.060257531 13 0.311582886 -0.467794233 -0.290852093 -1.443181647 14 -0.597049396 0.270528704 0.032229560 -1.013519905 15 3.750853466 0.037450201 -0.793536785 -0.510496279 16 -0.737115716 0.670544446 -1.104275776 -0.495770579 17 0.557994198 -0.035899739 0.451835868 0.037112527 18 0.340965301 0.653601911 0.516809543 1.904556951 19 1.434347302 1.909561116 -1.194707860 -0.313197177 > mylevel <- as.numeric(m$fitted[,'level']) > myslope <- as.numeric(m$fitted[,'slope']) > myseas <- as.numeric(m$fitted[,'sea']) > myresid <- as.numeric(m$resid) > myfit <- mylevel+myseas > mylagmax <- nx/2 > postscript(file="/var/www/html/freestat/rcomp/tmp/1flfu1291823130.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow = c(2,2)) > acf(as.numeric(x),lag.max = mylagmax,main='Observed') > acf(mylevel,na.action=na.pass,lag.max = mylagmax,main='Level') > acf(myseas,na.action=na.pass,lag.max = mylagmax,main='Seasonal') > acf(myresid,na.action=na.pass,lag.max = mylagmax,main='Standardized Residals') > par(op) > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/2pcwx1291823130.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow = c(2,2)) > spectrum(as.numeric(x),main='Observed') > spectrum(mylevel,main='Level') > spectrum(myseas,main='Seasonal') > spectrum(myresid,main='Standardized Residals') > par(op) > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/30mdi1291823130.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow = c(2,2)) > cpgram(as.numeric(x),main='Observed') > cpgram(mylevel,main='Level') > cpgram(myseas,main='Seasonal') > cpgram(myresid,main='Standardized Residals') > par(op) > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/40mdi1291823130.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > plot(as.numeric(m$resid),main='Standardized Residuals',ylab='Residuals',xlab='time',type='b') > grid() > dev.off() null device 1 > postscript(file="/var/www/html/freestat/rcomp/tmp/50mdi1291823130.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) > op <- par(mfrow = c(2,2)) > hist(m$resid,main='Residual Histogram') > plot(density(m$resid),main='Residual Kernel Density') > qqnorm(m$resid,main='Residual Normal QQ Plot') > qqline(m$resid) > plot(m$resid^2, myfit^2,main='Sq.Resid vs. Sq.Fit',xlab='Squared residuals',ylab='Squared Fit') > par(op) > dev.off() null device 1 > > #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab > load(file="/var/www/html/freestat/rcomp/createtable") > > a<-table.start() > a<-table.row.start(a) > a<-table.element(a,'Structural Time Series Model',6,TRUE) > a<-table.row.end(a) > a<-table.row.start(a) > a<-table.element(a,'t',header=TRUE) > a<-table.element(a,'Observed',header=TRUE) > a<-table.element(a,'Level',header=TRUE) > a<-table.element(a,'Slope',header=TRUE) > a<-table.element(a,'Seasonal',header=TRUE) > a<-table.element(a,'Stand. Residuals',header=TRUE) > a<-table.row.end(a) > for (i in 1:nx) { + a<-table.row.start(a) + a<-table.element(a,i,header=TRUE) + a<-table.element(a,x[i]) + a<-table.element(a,mylevel[i]) + a<-table.element(a,myslope[i]) + a<-table.element(a,myseas[i]) + a<-table.element(a,myresid[i]) + a<-table.row.end(a) + } > a<-table.end(a) > table.save(a,file="/var/www/html/freestat/rcomp/tmp/675su1291823130.tab") > > try(system("convert tmp/1flfu1291823130.ps tmp/1flfu1291823130.png",intern=TRUE)) character(0) > try(system("convert tmp/2pcwx1291823130.ps tmp/2pcwx1291823130.png",intern=TRUE)) character(0) > try(system("convert tmp/30mdi1291823130.ps tmp/30mdi1291823130.png",intern=TRUE)) character(0) > try(system("convert tmp/40mdi1291823130.ps tmp/40mdi1291823130.png",intern=TRUE)) character(0) > try(system("convert tmp/50mdi1291823130.ps tmp/50mdi1291823130.png",intern=TRUE)) character(0) > > > proc.time() user system elapsed 1.802 1.154 2.015