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(9782,9938,10111,10259,10419,10622,11173,11542,11538,11837,12060,12423,12791,12891,13098,13418,13614,13653,13980,14087,14332,14232,14226,14186,14310,14152,14127,14163,13964,13811,14440,14724,14790,14961,15117,15452,16080,16284,16524,16782,16663,16678,17448,17745,17789,17864,18079,18483,19037,19344,19590,19862,20207,20593,21253,21507,21528,21818,22205,22621,23006,23178,23358,23519,23725,23789,24472,24773,24477,24669,24827)
> par1 = '12'
> #'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 
8771.770 3695.355 1763.446    0.000 
> m$fitted
          level      slope         sea
Jan 1  9782.000   0.000000    0.000000
Feb 1  9918.936  23.554930   19.063942
Mar 1 10077.976  71.418823   33.023945
Apr 1 10241.984 110.157938   17.015584
May 1 10406.579 133.353645   12.421409
Jun 1 10604.508 161.004721   17.492018
Jul 1 11096.882 303.274224   76.118122
Aug 1 11530.253 359.188582   11.747307
Sep 1 11604.658 236.747128  -66.658468
Oct 1 11815.682 225.685804   21.317713
Nov 1 12049.720 229.277439   10.280180
Dec 1 12391.772 277.778087   31.228088
Jan 2 12802.435 333.404318  -11.434757
Feb 2 12943.231 254.498123  -52.231262
Mar 2 13101.989 215.731152   -3.988569
Apr 2 13383.943 243.427992   34.057344
May 2 13617.851 239.435723   -3.851152
Jun 2 13747.668 193.580505  -94.668077
Jul 2 13940.379 193.216765   39.620764
Aug 2 14044.956 156.102547   42.044243
Sep 2 14348.511 217.867216  -16.510588
Oct 2 14297.320 105.142821  -65.319504
Nov 2 14268.610  49.101914  -42.609911
Dec 2 14207.844   3.259785  -21.843993
Jan 3 14242.021  16.168527   67.979157
Feb 3 14196.361  -9.379154  -44.360955
Mar 3 14156.553 -21.771507  -29.553272
Apr 3 14127.583 -24.747002   35.417434
May 3 13975.505 -77.677647  -11.505223
Jun 3 13898.067 -77.578407  -87.066746
Jul 3 14264.558 106.584526  175.441747
Aug 3 14668.862 230.140402   55.138227
Sep 3 14793.267 186.232847   -3.266907
Oct 3 14990.608 190.845302  -29.607721
Nov 3 15147.531 176.778771  -30.531263
Dec 3 15446.666 227.425561    5.334060
Jan 4 15925.348 331.712868  154.651623
Feb 4 16314.396 355.372980  -30.395947
Mar 4 16570.784 314.913811  -46.784318
Apr 4 16734.566 252.738929   47.434117
May 4 16737.269 149.313078  -74.268901
Jun 4 16881.131 147.059879 -203.130934
Jul 4 17279.947 251.030099  168.053176
Aug 4 17639.521 295.873755  105.478855
Sep 4 17823.300 249.539611  -34.299797
Oct 4 17917.820 185.479044  -53.819967
Nov 4 18131.660 197.184896  -52.659816
Dec 4 18502.759 268.942952  -19.759478
Jan 5 18878.549 313.115456  158.450565
Feb 5 19317.519 364.949628   26.480648
Mar 5 19612.078 336.135117  -22.077753
Apr 5 19784.058 268.741773   77.942092
May 5 20230.122 341.856012  -23.121729
Jun 5 20814.753 441.969432 -221.752776
Jul 5 21155.657 400.331898   97.342895
Aug 5 21403.365 337.453653  103.634815
Sep 5 21558.744 262.422357  -30.743518
Oct 5 21859.057 278.032002  -41.057227
Nov 5 22268.239 332.020362  -63.238842
Dec 5 22659.767 356.526362  -38.766645
Jan 6 22908.815 312.211313   97.185215
Feb 6 23137.295 277.781569   40.705350
Mar 6 23350.014 251.137583    7.986225
Apr 6 23521.828 218.608602   -2.828088
May 6 23800.673 243.390825  -75.672581
Jun 6 23998.460 224.618264 -209.460288
Jul 6 24301.234 256.765218  170.765947
Aug 6 24606.473 276.696632  166.527298
Sep 6 24603.137 161.553761 -126.137276
Oct 6 24747.050 154.302213  -78.049784
Nov 6 24907.705 156.912893  -80.704987
> m$resid
           Jan          Feb          Mar          Apr          May          Jun
1  0.000000000  0.857758874  1.033918901  0.659212746  0.382204322  0.455522952
2  1.018222233 -1.219701814 -0.643610113  0.457873755 -0.065140367 -0.751501867
3  0.218397239 -0.412626620 -0.203793908 -0.049325742 -0.867776581  0.001626160
4  1.732823269  0.386059569 -0.664351703 -1.028947245 -1.700673363 -0.036951276
5  0.729674662  0.848644212 -0.473193731 -1.113284935  1.203915814  1.643496833
6 -0.730316512 -0.564538663 -0.437690459 -0.536673749  0.408283129 -0.308444981
           Jul          Aug          Sep          Oct          Nov          Dec
1  2.342187204  0.920043660 -2.014334518 -0.181964701  0.059083518  0.797847831
2 -0.005979392 -0.610490614  1.016044716 -1.854404585 -0.921850049 -0.756502856
3  3.025934577  2.032196301 -0.722291908  0.075875607 -0.231450508  0.836979559
4  1.707761458  0.737440087 -0.762167463 -1.053790378  0.192706238  1.185239195
5 -0.683867502 -1.033794483 -1.234104817  0.256800683  0.889103448  0.404409472
6  0.528045894  0.327635918 -1.893686720 -0.119314672  0.043001300             
> 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/13aai1293562635.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/2w1931293562635.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/3w1931293562635.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/4pa861293562635.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/5pa861293562635.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/63koe1293562635.tab") 
> 
> try(system("convert tmp/13aai1293562635.ps tmp/13aai1293562635.png",intern=TRUE))
character(0)
> try(system("convert tmp/2w1931293562635.ps tmp/2w1931293562635.png",intern=TRUE))
character(0)
> try(system("convert tmp/3w1931293562635.ps tmp/3w1931293562635.png",intern=TRUE))
character(0)
> try(system("convert tmp/4pa861293562635.ps tmp/4pa861293562635.png",intern=TRUE))
character(0)
> try(system("convert tmp/5pa861293562635.ps tmp/5pa861293562635.png",intern=TRUE))
character(0)
> 
> 
> proc.time()
   user  system elapsed 
  1.974   1.169   2.140