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.
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 <- array(list(32.68
+ ,10967.87
+ ,31.54
+ ,10433.56
+ ,32.43
+ ,10665.78
+ ,26.54
+ ,10666.71
+ ,25.85
+ ,10682.74
+ ,27.6
+ ,10777.22
+ ,25.71
+ ,10052.6
+ ,25.38
+ ,10213.97
+ ,28.57
+ ,10546.82
+ ,27.64
+ ,10767.2
+ ,25.36
+ ,10444.5
+ ,25.9
+ ,10314.68
+ ,26.29
+ ,9042.56
+ ,21.74
+ ,9220.75
+ ,19.2
+ ,9721.84
+ ,19.32
+ ,9978.53
+ ,19.82
+ ,9923.81
+ ,20.36
+ ,9892.56
+ ,24.31
+ ,10500.98
+ ,25.97
+ ,10179.35
+ ,25.61
+ ,10080.48
+ ,24.67
+ ,9492.44
+ ,25.59
+ ,8616.49
+ ,26.09
+ ,8685.4
+ ,28.37
+ ,8160.67
+ ,27.34
+ ,8048.1
+ ,24.46
+ ,8641.21
+ ,27.46
+ ,8526.63
+ ,30.23
+ ,8474.21
+ ,32.33
+ ,7916.13
+ ,29.87
+ ,7977.64
+ ,24.87
+ ,8334.59
+ ,25.48
+ ,8623.36
+ ,27.28
+ ,9098.03
+ ,28.24
+ ,9154.34
+ ,29.58
+ ,9284.73
+ ,26.95
+ ,9492.49
+ ,29.08
+ ,9682.35
+ ,28.76
+ ,9762.12
+ ,29.59
+ ,10124.63
+ ,30.7
+ ,10540.05
+ ,30.52
+ ,10601.61
+ ,32.67
+ ,10323.73
+ ,33.19
+ ,10418.4
+ ,37.13
+ ,10092.96
+ ,35.54
+ ,10364.91
+ ,37.75
+ ,10152.09
+ ,41.84
+ ,10032.8
+ ,42.94
+ ,10204.59
+ ,49.14
+ ,10001.6
+ ,44.61
+ ,10411.75
+ ,40.22
+ ,10673.38
+ ,44.23
+ ,10539.51
+ ,45.85
+ ,10723.78
+ ,53.38
+ ,10682.06
+ ,53.26
+ ,10283.19
+ ,51.8
+ ,10377.18
+ ,55.3
+ ,10486.64
+ ,57.81
+ ,10545.38
+ ,63.96
+ ,10554.27
+ ,63.77
+ ,10532.54
+ ,59.15
+ ,10324.31
+ ,56.12
+ ,10695.25
+ ,57.42
+ ,10827.81
+ ,63.52
+ ,10872.48
+ ,61.71
+ ,10971.19
+ ,63.01
+ ,11145.65
+ ,68.18
+ ,11234.68
+ ,72.03
+ ,11333.88
+ ,69.75
+ ,10997.97
+ ,74.41
+ ,11036.89
+ ,74.33
+ ,11257.35
+ ,64.24
+ ,11533.59
+ ,60.03
+ ,11963.12
+ ,59.44
+ ,12185.15
+ ,62.5
+ ,12377.62
+ ,55.04
+ ,12512.89
+ ,58.34
+ ,12631.48
+ ,61.92
+ ,12268.53
+ ,67.65
+ ,12754.8
+ ,67.68
+ ,13407.75
+ ,70.3
+ ,13480.21
+ ,75.26
+ ,13673.28
+ ,71.44
+ ,13239.71
+ ,76.36
+ ,13557.69
+ ,81.71
+ ,13901.28
+ ,92.6
+ ,13200.58
+ ,90.6
+ ,13406.97
+ ,92.23
+ ,12538.12
+ ,94.09
+ ,12419.57
+ ,102.79
+ ,12193.88
+ ,109.65
+ ,12656.63
+ ,124.05
+ ,12812.48
+ ,132.69
+ ,12056.67
+ ,135.81
+ ,11322.38
+ ,116.07
+ ,11530.75
+ ,101.42
+ ,11114.08
+ ,75.73
+ ,9181.73
+ ,55.48
+ ,8614.55)
+ ,dim=c(2
+ ,99)
+ ,dimnames=list(c('Olieprijs'
+ ,'DowJones')
+ ,1:99))
> y <- array(NA,dim=c(2,99),dimnames=list(c('Olieprijs','DowJones'),1:99))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'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!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Olieprijs DowJones M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 32.68 10967.87 1 0 0 0 0 0 0 0 0 0 0
2 31.54 10433.56 0 1 0 0 0 0 0 0 0 0 0
3 32.43 10665.78 0 0 1 0 0 0 0 0 0 0 0
4 26.54 10666.71 0 0 0 1 0 0 0 0 0 0 0
5 25.85 10682.74 0 0 0 0 1 0 0 0 0 0 0
6 27.60 10777.22 0 0 0 0 0 1 0 0 0 0 0
7 25.71 10052.60 0 0 0 0 0 0 1 0 0 0 0
8 25.38 10213.97 0 0 0 0 0 0 0 1 0 0 0
9 28.57 10546.82 0 0 0 0 0 0 0 0 1 0 0
10 27.64 10767.20 0 0 0 0 0 0 0 0 0 1 0
11 25.36 10444.50 0 0 0 0 0 0 0 0 0 0 1
12 25.90 10314.68 0 0 0 0 0 0 0 0 0 0 0
13 26.29 9042.56 1 0 0 0 0 0 0 0 0 0 0
14 21.74 9220.75 0 1 0 0 0 0 0 0 0 0 0
15 19.20 9721.84 0 0 1 0 0 0 0 0 0 0 0
16 19.32 9978.53 0 0 0 1 0 0 0 0 0 0 0
17 19.82 9923.81 0 0 0 0 1 0 0 0 0 0 0
18 20.36 9892.56 0 0 0 0 0 1 0 0 0 0 0
19 24.31 10500.98 0 0 0 0 0 0 1 0 0 0 0
20 25.97 10179.35 0 0 0 0 0 0 0 1 0 0 0
21 25.61 10080.48 0 0 0 0 0 0 0 0 1 0 0
22 24.67 9492.44 0 0 0 0 0 0 0 0 0 1 0
23 25.59 8616.49 0 0 0 0 0 0 0 0 0 0 1
24 26.09 8685.40 0 0 0 0 0 0 0 0 0 0 0
25 28.37 8160.67 1 0 0 0 0 0 0 0 0 0 0
26 27.34 8048.10 0 1 0 0 0 0 0 0 0 0 0
27 24.46 8641.21 0 0 1 0 0 0 0 0 0 0 0
28 27.46 8526.63 0 0 0 1 0 0 0 0 0 0 0
29 30.23 8474.21 0 0 0 0 1 0 0 0 0 0 0
30 32.33 7916.13 0 0 0 0 0 1 0 0 0 0 0
31 29.87 7977.64 0 0 0 0 0 0 1 0 0 0 0
32 24.87 8334.59 0 0 0 0 0 0 0 1 0 0 0
33 25.48 8623.36 0 0 0 0 0 0 0 0 1 0 0
34 27.28 9098.03 0 0 0 0 0 0 0 0 0 1 0
35 28.24 9154.34 0 0 0 0 0 0 0 0 0 0 1
36 29.58 9284.73 0 0 0 0 0 0 0 0 0 0 0
37 26.95 9492.49 1 0 0 0 0 0 0 0 0 0 0
38 29.08 9682.35 0 1 0 0 0 0 0 0 0 0 0
39 28.76 9762.12 0 0 1 0 0 0 0 0 0 0 0
40 29.59 10124.63 0 0 0 1 0 0 0 0 0 0 0
41 30.70 10540.05 0 0 0 0 1 0 0 0 0 0 0
42 30.52 10601.61 0 0 0 0 0 1 0 0 0 0 0
43 32.67 10323.73 0 0 0 0 0 0 1 0 0 0 0
44 33.19 10418.40 0 0 0 0 0 0 0 1 0 0 0
45 37.13 10092.96 0 0 0 0 0 0 0 0 1 0 0
46 35.54 10364.91 0 0 0 0 0 0 0 0 0 1 0
47 37.75 10152.09 0 0 0 0 0 0 0 0 0 0 1
48 41.84 10032.80 0 0 0 0 0 0 0 0 0 0 0
49 42.94 10204.59 1 0 0 0 0 0 0 0 0 0 0
50 49.14 10001.60 0 1 0 0 0 0 0 0 0 0 0
51 44.61 10411.75 0 0 1 0 0 0 0 0 0 0 0
52 40.22 10673.38 0 0 0 1 0 0 0 0 0 0 0
53 44.23 10539.51 0 0 0 0 1 0 0 0 0 0 0
54 45.85 10723.78 0 0 0 0 0 1 0 0 0 0 0
55 53.38 10682.06 0 0 0 0 0 0 1 0 0 0 0
56 53.26 10283.19 0 0 0 0 0 0 0 1 0 0 0
57 51.80 10377.18 0 0 0 0 0 0 0 0 1 0 0
58 55.30 10486.64 0 0 0 0 0 0 0 0 0 1 0
59 57.81 10545.38 0 0 0 0 0 0 0 0 0 0 1
60 63.96 10554.27 0 0 0 0 0 0 0 0 0 0 0
61 63.77 10532.54 1 0 0 0 0 0 0 0 0 0 0
62 59.15 10324.31 0 1 0 0 0 0 0 0 0 0 0
63 56.12 10695.25 0 0 1 0 0 0 0 0 0 0 0
64 57.42 10827.81 0 0 0 1 0 0 0 0 0 0 0
65 63.52 10872.48 0 0 0 0 1 0 0 0 0 0 0
66 61.71 10971.19 0 0 0 0 0 1 0 0 0 0 0
67 63.01 11145.65 0 0 0 0 0 0 1 0 0 0 0
68 68.18 11234.68 0 0 0 0 0 0 0 1 0 0 0
69 72.03 11333.88 0 0 0 0 0 0 0 0 1 0 0
70 69.75 10997.97 0 0 0 0 0 0 0 0 0 1 0
71 74.41 11036.89 0 0 0 0 0 0 0 0 0 0 1
72 74.33 11257.35 0 0 0 0 0 0 0 0 0 0 0
73 64.24 11533.59 1 0 0 0 0 0 0 0 0 0 0
74 60.03 11963.12 0 1 0 0 0 0 0 0 0 0 0
75 59.44 12185.15 0 0 1 0 0 0 0 0 0 0 0
76 62.50 12377.62 0 0 0 1 0 0 0 0 0 0 0
77 55.04 12512.89 0 0 0 0 1 0 0 0 0 0 0
78 58.34 12631.48 0 0 0 0 0 1 0 0 0 0 0
79 61.92 12268.53 0 0 0 0 0 0 1 0 0 0 0
80 67.65 12754.80 0 0 0 0 0 0 0 1 0 0 0
81 67.68 13407.75 0 0 0 0 0 0 0 0 1 0 0
82 70.30 13480.21 0 0 0 0 0 0 0 0 0 1 0
83 75.26 13673.28 0 0 0 0 0 0 0 0 0 0 1
84 71.44 13239.71 0 0 0 0 0 0 0 0 0 0 0
85 76.36 13557.69 1 0 0 0 0 0 0 0 0 0 0
86 81.71 13901.28 0 1 0 0 0 0 0 0 0 0 0
87 92.60 13200.58 0 0 1 0 0 0 0 0 0 0 0
88 90.60 13406.97 0 0 0 1 0 0 0 0 0 0 0
89 92.23 12538.12 0 0 0 0 1 0 0 0 0 0 0
90 94.09 12419.57 0 0 0 0 0 1 0 0 0 0 0
91 102.79 12193.88 0 0 0 0 0 0 1 0 0 0 0
92 109.65 12656.63 0 0 0 0 0 0 0 1 0 0 0
93 124.05 12812.48 0 0 0 0 0 0 0 0 1 0 0
94 132.69 12056.67 0 0 0 0 0 0 0 0 0 1 0
95 135.81 11322.38 0 0 0 0 0 0 0 0 0 0 1
96 116.07 11530.75 0 0 0 0 0 0 0 0 0 0 0
97 101.42 11114.08 1 0 0 0 0 0 0 0 0 0 0
98 75.73 9181.73 0 1 0 0 0 0 0 0 0 0 0
99 55.48 8614.55 0 0 1 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) DowJones M1 M2 M3 M4
-86.28428 0.01342 -3.35338 -3.65780 -7.84444 -14.76787
M5 M6 M7 M8 M9 M10
-12.93534 -11.53573 -7.35542 -7.10533 -6.09239 -3.84930
M11
1.30090
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-29.838 -13.816 -4.696 10.040 68.830
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -86.284276 17.401376 -4.958 3.53e-06 ***
DowJones 0.013422 0.001488 9.022 4.45e-14 ***
M1 -3.353376 10.057783 -0.333 0.740
M2 -3.657801 10.066976 -0.363 0.717
M3 -7.844439 10.060205 -0.780 0.438
M4 -14.767868 10.352955 -1.426 0.157
M5 -12.935339 10.350568 -1.250 0.215
M6 -11.535732 10.350011 -1.115 0.268
M7 -7.355422 10.348326 -0.711 0.479
M8 -7.105332 10.350536 -0.686 0.494
M9 -6.092387 10.357648 -0.588 0.558
M10 -3.849301 10.353908 -0.372 0.711
M11 1.300897 10.348229 0.126 0.900
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.7 on 86 degrees of freedom
Multiple R-squared: 0.4993, Adjusted R-squared: 0.4295
F-statistic: 7.147 on 12 and 86 DF, p-value: 7.205e-09
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 1.041616e-02 2.083232e-02 0.9895838
[2,] 1.735292e-03 3.470585e-03 0.9982647
[3,] 2.717790e-04 5.435579e-04 0.9997282
[4,] 6.481139e-05 1.296228e-04 0.9999352
[5,] 9.594350e-06 1.918870e-05 0.9999904
[6,] 1.334213e-06 2.668425e-06 0.9999987
[7,] 3.695137e-07 7.390275e-07 0.9999996
[8,] 4.771146e-07 9.542292e-07 0.9999995
[9,] 1.844862e-07 3.689724e-07 0.9999998
[10,] 5.390698e-08 1.078140e-07 0.9999999
[11,] 1.793426e-08 3.586852e-08 1.0000000
[12,] 3.277343e-09 6.554686e-09 1.0000000
[13,] 2.401862e-09 4.803724e-09 1.0000000
[14,] 2.802566e-09 5.605131e-09 1.0000000
[15,] 2.747742e-09 5.495483e-09 1.0000000
[16,] 7.979371e-10 1.595874e-09 1.0000000
[17,] 1.568120e-10 3.136239e-10 1.0000000
[18,] 3.137138e-11 6.274275e-11 1.0000000
[19,] 6.247310e-12 1.249462e-11 1.0000000
[20,] 1.481424e-12 2.962847e-12 1.0000000
[21,] 3.810222e-13 7.620444e-13 1.0000000
[22,] 8.576516e-14 1.715303e-13 1.0000000
[23,] 1.929550e-14 3.859100e-14 1.0000000
[24,] 4.948038e-15 9.896075e-15 1.0000000
[25,] 1.721716e-15 3.443433e-15 1.0000000
[26,] 7.071510e-16 1.414302e-15 1.0000000
[27,] 2.078000e-16 4.156001e-16 1.0000000
[28,] 1.066804e-16 2.133608e-16 1.0000000
[29,] 1.001733e-16 2.003465e-16 1.0000000
[30,] 2.308612e-16 4.617225e-16 1.0000000
[31,] 3.534186e-16 7.068371e-16 1.0000000
[32,] 1.162587e-15 2.325173e-15 1.0000000
[33,] 7.744863e-15 1.548973e-14 1.0000000
[34,] 2.972627e-14 5.945255e-14 1.0000000
[35,] 8.772544e-13 1.754509e-12 1.0000000
[36,] 3.730308e-12 7.460616e-12 1.0000000
[37,] 4.990301e-12 9.980602e-12 1.0000000
[38,] 1.156193e-11 2.312386e-11 1.0000000
[39,] 2.226500e-11 4.453001e-11 1.0000000
[40,] 1.447698e-10 2.895396e-10 1.0000000
[41,] 9.758768e-10 1.951754e-09 1.0000000
[42,] 3.140195e-09 6.280391e-09 1.0000000
[43,] 1.730283e-08 3.460566e-08 1.0000000
[44,] 9.234035e-08 1.846807e-07 0.9999999
[45,] 4.031280e-07 8.062560e-07 0.9999996
[46,] 1.065583e-06 2.131166e-06 0.9999989
[47,] 1.410712e-06 2.821423e-06 0.9999986
[48,] 1.580866e-06 3.161732e-06 0.9999984
[49,] 2.140932e-06 4.281864e-06 0.9999979
[50,] 3.373878e-06 6.747755e-06 0.9999966
[51,] 3.695168e-06 7.390337e-06 0.9999963
[52,] 3.643359e-06 7.286718e-06 0.9999964
[53,] 4.651529e-06 9.303059e-06 0.9999953
[54,] 7.003362e-06 1.400672e-05 0.9999930
[55,] 1.506792e-05 3.013583e-05 0.9999849
[56,] 3.245208e-05 6.490417e-05 0.9999675
[57,] 3.267410e-05 6.534820e-05 0.9999673
[58,] 2.676189e-05 5.352378e-05 0.9999732
[59,] 1.452080e-05 2.904159e-05 0.9999855
[60,] 6.528752e-06 1.305750e-05 0.9999935
[61,] 4.920679e-06 9.841358e-06 0.9999951
[62,] 4.918299e-06 9.836599e-06 0.9999951
[63,] 4.395701e-06 8.791402e-06 0.9999956
[64,] 5.401592e-06 1.080318e-05 0.9999946
[65,] 7.520335e-06 1.504067e-05 0.9999925
[66,] 3.618901e-05 7.237801e-05 0.9999638
[67,] 4.120394e-04 8.240788e-04 0.9995880
[68,] 9.409236e-03 1.881847e-02 0.9905908
> postscript(file="/var/www/html/rcomp/tmp/19ag81229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2nwki1229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3x3f71229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4ljdd1229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5i8n61229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 99
Frequency = 1
1 2 3 4 5 6
-24.8880053 -18.5523199 -16.5924308 -15.5714838 -18.3091602 -19.2268336
7 8 9 10 11 12
-15.5716303 -18.3175540 -20.6078562 -26.7387802 -29.8378487 -26.2545684
13 14 15 16 17 18
-5.4373888 -12.0745474 -17.1533065 -13.5550515 -14.1531542 -14.5933380
19 20 21 22 23 24
-22.9895784 -17.2629005 -17.3088573 -12.5995435 -5.0731475 -4.1971290
25 26 27 28 29 30
8.4789292 9.2642162 2.6104074 14.0716759 15.7127037 23.9033875
31 32 33 34 35 36
16.4375191 6.3966114 2.1179306 -4.6959557 -9.6419204 -8.7510579
37 38 39 40 41 42
-10.8161402 -10.9299282 -8.1339260 -5.2459379 -11.5440413 -13.9498777
43 44 45 46 47 48
-12.2506112 -13.2513185 -5.9563581 -13.4394306 -13.5232575 -6.5313059
49 50 51 52 53 54
-4.3836157 4.8452466 -1.0029581 -1.9810055 1.9932064 -0.2595867
55 56 57 58 59 60
3.6500500 8.6334073 4.8989731 4.6867658 1.2581868 8.5897659
61 62 63 64 65 66
12.0447917 10.5239828 6.7020366 13.1463068 16.8142377 12.2797911
67 68 69 70 71 72
7.0579603 10.7829504 12.2885901 12.2739321 11.2613679 9.5233527
73 74 75 76 77 78
-0.9208364 -10.5913638 -9.9747092 -2.5745232 -13.6825834 -13.3738501
79 80 81 82 83 84
-9.1028135 -10.1493941 -19.8959298 -20.4915403 -23.2730342 -19.9729624
85 86 87 88 89 90
-15.9673665 -14.9244471 9.5566610 11.7100192 23.1687913 25.2203075
91 92 93 94 95 96
32.7691041 33.1681980 44.4635076 61.0045524 68.8296537 47.5939049
97 98 99
41.8896320 42.4391609 33.9882256
> postscript(file="/var/www/html/rcomp/tmp/6loo31229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 99
Frequency = 1
lag(myerror, k = 1) myerror
0 -24.8880053 NA
1 -18.5523199 -24.8880053
2 -16.5924308 -18.5523199
3 -15.5714838 -16.5924308
4 -18.3091602 -15.5714838
5 -19.2268336 -18.3091602
6 -15.5716303 -19.2268336
7 -18.3175540 -15.5716303
8 -20.6078562 -18.3175540
9 -26.7387802 -20.6078562
10 -29.8378487 -26.7387802
11 -26.2545684 -29.8378487
12 -5.4373888 -26.2545684
13 -12.0745474 -5.4373888
14 -17.1533065 -12.0745474
15 -13.5550515 -17.1533065
16 -14.1531542 -13.5550515
17 -14.5933380 -14.1531542
18 -22.9895784 -14.5933380
19 -17.2629005 -22.9895784
20 -17.3088573 -17.2629005
21 -12.5995435 -17.3088573
22 -5.0731475 -12.5995435
23 -4.1971290 -5.0731475
24 8.4789292 -4.1971290
25 9.2642162 8.4789292
26 2.6104074 9.2642162
27 14.0716759 2.6104074
28 15.7127037 14.0716759
29 23.9033875 15.7127037
30 16.4375191 23.9033875
31 6.3966114 16.4375191
32 2.1179306 6.3966114
33 -4.6959557 2.1179306
34 -9.6419204 -4.6959557
35 -8.7510579 -9.6419204
36 -10.8161402 -8.7510579
37 -10.9299282 -10.8161402
38 -8.1339260 -10.9299282
39 -5.2459379 -8.1339260
40 -11.5440413 -5.2459379
41 -13.9498777 -11.5440413
42 -12.2506112 -13.9498777
43 -13.2513185 -12.2506112
44 -5.9563581 -13.2513185
45 -13.4394306 -5.9563581
46 -13.5232575 -13.4394306
47 -6.5313059 -13.5232575
48 -4.3836157 -6.5313059
49 4.8452466 -4.3836157
50 -1.0029581 4.8452466
51 -1.9810055 -1.0029581
52 1.9932064 -1.9810055
53 -0.2595867 1.9932064
54 3.6500500 -0.2595867
55 8.6334073 3.6500500
56 4.8989731 8.6334073
57 4.6867658 4.8989731
58 1.2581868 4.6867658
59 8.5897659 1.2581868
60 12.0447917 8.5897659
61 10.5239828 12.0447917
62 6.7020366 10.5239828
63 13.1463068 6.7020366
64 16.8142377 13.1463068
65 12.2797911 16.8142377
66 7.0579603 12.2797911
67 10.7829504 7.0579603
68 12.2885901 10.7829504
69 12.2739321 12.2885901
70 11.2613679 12.2739321
71 9.5233527 11.2613679
72 -0.9208364 9.5233527
73 -10.5913638 -0.9208364
74 -9.9747092 -10.5913638
75 -2.5745232 -9.9747092
76 -13.6825834 -2.5745232
77 -13.3738501 -13.6825834
78 -9.1028135 -13.3738501
79 -10.1493941 -9.1028135
80 -19.8959298 -10.1493941
81 -20.4915403 -19.8959298
82 -23.2730342 -20.4915403
83 -19.9729624 -23.2730342
84 -15.9673665 -19.9729624
85 -14.9244471 -15.9673665
86 9.5566610 -14.9244471
87 11.7100192 9.5566610
88 23.1687913 11.7100192
89 25.2203075 23.1687913
90 32.7691041 25.2203075
91 33.1681980 32.7691041
92 44.4635076 33.1681980
93 61.0045524 44.4635076
94 68.8296537 61.0045524
95 47.5939049 68.8296537
96 41.8896320 47.5939049
97 42.4391609 41.8896320
98 33.9882256 42.4391609
99 NA 33.9882256
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -18.5523199 -24.8880053
[2,] -16.5924308 -18.5523199
[3,] -15.5714838 -16.5924308
[4,] -18.3091602 -15.5714838
[5,] -19.2268336 -18.3091602
[6,] -15.5716303 -19.2268336
[7,] -18.3175540 -15.5716303
[8,] -20.6078562 -18.3175540
[9,] -26.7387802 -20.6078562
[10,] -29.8378487 -26.7387802
[11,] -26.2545684 -29.8378487
[12,] -5.4373888 -26.2545684
[13,] -12.0745474 -5.4373888
[14,] -17.1533065 -12.0745474
[15,] -13.5550515 -17.1533065
[16,] -14.1531542 -13.5550515
[17,] -14.5933380 -14.1531542
[18,] -22.9895784 -14.5933380
[19,] -17.2629005 -22.9895784
[20,] -17.3088573 -17.2629005
[21,] -12.5995435 -17.3088573
[22,] -5.0731475 -12.5995435
[23,] -4.1971290 -5.0731475
[24,] 8.4789292 -4.1971290
[25,] 9.2642162 8.4789292
[26,] 2.6104074 9.2642162
[27,] 14.0716759 2.6104074
[28,] 15.7127037 14.0716759
[29,] 23.9033875 15.7127037
[30,] 16.4375191 23.9033875
[31,] 6.3966114 16.4375191
[32,] 2.1179306 6.3966114
[33,] -4.6959557 2.1179306
[34,] -9.6419204 -4.6959557
[35,] -8.7510579 -9.6419204
[36,] -10.8161402 -8.7510579
[37,] -10.9299282 -10.8161402
[38,] -8.1339260 -10.9299282
[39,] -5.2459379 -8.1339260
[40,] -11.5440413 -5.2459379
[41,] -13.9498777 -11.5440413
[42,] -12.2506112 -13.9498777
[43,] -13.2513185 -12.2506112
[44,] -5.9563581 -13.2513185
[45,] -13.4394306 -5.9563581
[46,] -13.5232575 -13.4394306
[47,] -6.5313059 -13.5232575
[48,] -4.3836157 -6.5313059
[49,] 4.8452466 -4.3836157
[50,] -1.0029581 4.8452466
[51,] -1.9810055 -1.0029581
[52,] 1.9932064 -1.9810055
[53,] -0.2595867 1.9932064
[54,] 3.6500500 -0.2595867
[55,] 8.6334073 3.6500500
[56,] 4.8989731 8.6334073
[57,] 4.6867658 4.8989731
[58,] 1.2581868 4.6867658
[59,] 8.5897659 1.2581868
[60,] 12.0447917 8.5897659
[61,] 10.5239828 12.0447917
[62,] 6.7020366 10.5239828
[63,] 13.1463068 6.7020366
[64,] 16.8142377 13.1463068
[65,] 12.2797911 16.8142377
[66,] 7.0579603 12.2797911
[67,] 10.7829504 7.0579603
[68,] 12.2885901 10.7829504
[69,] 12.2739321 12.2885901
[70,] 11.2613679 12.2739321
[71,] 9.5233527 11.2613679
[72,] -0.9208364 9.5233527
[73,] -10.5913638 -0.9208364
[74,] -9.9747092 -10.5913638
[75,] -2.5745232 -9.9747092
[76,] -13.6825834 -2.5745232
[77,] -13.3738501 -13.6825834
[78,] -9.1028135 -13.3738501
[79,] -10.1493941 -9.1028135
[80,] -19.8959298 -10.1493941
[81,] -20.4915403 -19.8959298
[82,] -23.2730342 -20.4915403
[83,] -19.9729624 -23.2730342
[84,] -15.9673665 -19.9729624
[85,] -14.9244471 -15.9673665
[86,] 9.5566610 -14.9244471
[87,] 11.7100192 9.5566610
[88,] 23.1687913 11.7100192
[89,] 25.2203075 23.1687913
[90,] 32.7691041 25.2203075
[91,] 33.1681980 32.7691041
[92,] 44.4635076 33.1681980
[93,] 61.0045524 44.4635076
[94,] 68.8296537 61.0045524
[95,] 47.5939049 68.8296537
[96,] 41.8896320 47.5939049
[97,] 42.4391609 41.8896320
[98,] 33.9882256 42.4391609
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -18.5523199 -24.8880053
2 -16.5924308 -18.5523199
3 -15.5714838 -16.5924308
4 -18.3091602 -15.5714838
5 -19.2268336 -18.3091602
6 -15.5716303 -19.2268336
7 -18.3175540 -15.5716303
8 -20.6078562 -18.3175540
9 -26.7387802 -20.6078562
10 -29.8378487 -26.7387802
11 -26.2545684 -29.8378487
12 -5.4373888 -26.2545684
13 -12.0745474 -5.4373888
14 -17.1533065 -12.0745474
15 -13.5550515 -17.1533065
16 -14.1531542 -13.5550515
17 -14.5933380 -14.1531542
18 -22.9895784 -14.5933380
19 -17.2629005 -22.9895784
20 -17.3088573 -17.2629005
21 -12.5995435 -17.3088573
22 -5.0731475 -12.5995435
23 -4.1971290 -5.0731475
24 8.4789292 -4.1971290
25 9.2642162 8.4789292
26 2.6104074 9.2642162
27 14.0716759 2.6104074
28 15.7127037 14.0716759
29 23.9033875 15.7127037
30 16.4375191 23.9033875
31 6.3966114 16.4375191
32 2.1179306 6.3966114
33 -4.6959557 2.1179306
34 -9.6419204 -4.6959557
35 -8.7510579 -9.6419204
36 -10.8161402 -8.7510579
37 -10.9299282 -10.8161402
38 -8.1339260 -10.9299282
39 -5.2459379 -8.1339260
40 -11.5440413 -5.2459379
41 -13.9498777 -11.5440413
42 -12.2506112 -13.9498777
43 -13.2513185 -12.2506112
44 -5.9563581 -13.2513185
45 -13.4394306 -5.9563581
46 -13.5232575 -13.4394306
47 -6.5313059 -13.5232575
48 -4.3836157 -6.5313059
49 4.8452466 -4.3836157
50 -1.0029581 4.8452466
51 -1.9810055 -1.0029581
52 1.9932064 -1.9810055
53 -0.2595867 1.9932064
54 3.6500500 -0.2595867
55 8.6334073 3.6500500
56 4.8989731 8.6334073
57 4.6867658 4.8989731
58 1.2581868 4.6867658
59 8.5897659 1.2581868
60 12.0447917 8.5897659
61 10.5239828 12.0447917
62 6.7020366 10.5239828
63 13.1463068 6.7020366
64 16.8142377 13.1463068
65 12.2797911 16.8142377
66 7.0579603 12.2797911
67 10.7829504 7.0579603
68 12.2885901 10.7829504
69 12.2739321 12.2885901
70 11.2613679 12.2739321
71 9.5233527 11.2613679
72 -0.9208364 9.5233527
73 -10.5913638 -0.9208364
74 -9.9747092 -10.5913638
75 -2.5745232 -9.9747092
76 -13.6825834 -2.5745232
77 -13.3738501 -13.6825834
78 -9.1028135 -13.3738501
79 -10.1493941 -9.1028135
80 -19.8959298 -10.1493941
81 -20.4915403 -19.8959298
82 -23.2730342 -20.4915403
83 -19.9729624 -23.2730342
84 -15.9673665 -19.9729624
85 -14.9244471 -15.9673665
86 9.5566610 -14.9244471
87 11.7100192 9.5566610
88 23.1687913 11.7100192
89 25.2203075 23.1687913
90 32.7691041 25.2203075
91 33.1681980 32.7691041
92 44.4635076 33.1681980
93 61.0045524 44.4635076
94 68.8296537 61.0045524
95 47.5939049 68.8296537
96 41.8896320 47.5939049
97 42.4391609 41.8896320
98 33.9882256 42.4391609
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/77an71229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8whb41229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9vfdr1229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10to4j1229181967.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/111dqm1229181967.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/120ye01229181967.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13k1x01229181967.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14cj7c1229181967.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15u1x91229181967.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/161fk11229181967.tab")
+ }
>
> system("convert tmp/19ag81229181967.ps tmp/19ag81229181967.png")
> system("convert tmp/2nwki1229181967.ps tmp/2nwki1229181967.png")
> system("convert tmp/3x3f71229181967.ps tmp/3x3f71229181967.png")
> system("convert tmp/4ljdd1229181967.ps tmp/4ljdd1229181967.png")
> system("convert tmp/5i8n61229181967.ps tmp/5i8n61229181967.png")
> system("convert tmp/6loo31229181967.ps tmp/6loo31229181967.png")
> system("convert tmp/77an71229181967.ps tmp/77an71229181967.png")
> system("convert tmp/8whb41229181967.ps tmp/8whb41229181967.png")
> system("convert tmp/9vfdr1229181967.ps tmp/9vfdr1229181967.png")
> system("convert tmp/10to4j1229181967.ps tmp/10to4j1229181967.png")
>
>
> proc.time()
user system elapsed
2.981 1.612 3.587