R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(10284.5
+ ,1.038351422
+ ,1.4
+ ,12792
+ ,0.933031106
+ ,1.3
+ ,12823.61538
+ ,0.932783124
+ ,1.3
+ ,13845.66667
+ ,0.953755367
+ ,1.2
+ ,15335.63636
+ ,1.009865664
+ ,1.1
+ ,11188.5
+ ,0.979532493
+ ,1.4
+ ,13633.25
+ ,0.98651077
+ ,1.2
+ ,12298.46667
+ ,0.964661281
+ ,1.5
+ ,15353.63636
+ ,0.946761816
+ ,1.1
+ ,12696.15385
+ ,0.959068881
+ ,1.3
+ ,12213.93333
+ ,0.985710058
+ ,1.5
+ ,13683.72727
+ ,0.92582159
+ ,1.1
+ ,11214.14286
+ ,1.036865325
+ ,1.4
+ ,13950.23077
+ ,0.944443576
+ ,1.3
+ ,11179.13333
+ ,0.944901812
+ ,1.5
+ ,11801.875
+ ,0.989151445
+ ,1.6
+ ,11188.82353
+ ,1.054361624
+ ,1.7
+ ,16456.27273
+ ,1.033478919
+ ,1.1
+ ,11110.0625
+ ,1.001368875
+ ,1.6
+ ,16530.69231
+ ,1.019812646
+ ,1.3
+ ,10038.41176
+ ,0.993902155
+ ,1.7
+ ,11681.25
+ ,0.961444482
+ ,1.6
+ ,11148.88235
+ ,0.957449711
+ ,1.7
+ ,8631
+ ,0.93308639
+ ,1.9
+ ,9386.444444
+ ,1.045170549
+ ,1.8
+ ,9764.736842
+ ,0.953166261
+ ,1.9
+ ,12043.75
+ ,0.966782226
+ ,1.6
+ ,12948.06667
+ ,0.972992606
+ ,1.5
+ ,10987.125
+ ,1.013607482
+ ,1.6
+ ,11648.3125
+ ,0.984839518
+ ,1.6
+ ,10633.35294
+ ,0.973220775
+ ,1.7
+ ,10219.3
+ ,0.957284573
+ ,2
+ ,9037.6
+ ,0.972067159
+ ,2
+ ,10296.31579
+ ,0.986878944
+ ,1.9
+ ,11705.41176
+ ,0.954654488
+ ,1.7
+ ,10681.94444
+ ,0.978986976
+ ,1.8
+ ,9362.947368
+ ,1.003056035
+ ,1.9
+ ,11306.35294
+ ,0.970081156
+ ,1.7
+ ,10984.45
+ ,0.991376354
+ ,2
+ ,10062.61905
+ ,1.022609041
+ ,2.1
+ ,8118.583333
+ ,1.089901216
+ ,2.4
+ ,8867.48
+ ,1.060373568
+ ,2.5
+ ,8346.72
+ ,0.985952627
+ ,2.5
+ ,8529.307692
+ ,1.037512164
+ ,2.6
+ ,10697.18182
+ ,1.025335152
+ ,2.2
+ ,8591.84
+ ,1.006376649
+ ,2.5
+ ,8695.607143
+ ,1.018762056
+ ,2.8
+ ,8125.571429
+ ,1.01601847
+ ,2.8
+ ,7009.758621
+ ,1.112410461
+ ,2.9
+ ,7883.466667
+ ,1.037903689
+ ,3
+ ,7527.645161
+ ,1.045436015
+ ,3.1
+ ,6763.758621
+ ,1.09935434
+ ,2.9
+ ,6682.333333
+ ,1.101920787
+ ,2.7
+ ,7855.681818
+ ,1.080574973
+ ,2.2
+ ,6738.88
+ ,1.024388761
+ ,2.5
+ ,7895.434783
+ ,1.024282249
+ ,2.3
+ ,6361.884615
+ ,0.993865289
+ ,2.6
+ ,6935.956522
+ ,0.984935203
+ ,2.3
+ ,8344.454545
+ ,1.005791114
+ ,2.2
+ ,9107.944444
+ ,0.94742834
+ ,1.8)
+ ,dim=c(3
+ ,60)
+ ,dimnames=list(c('Invoer/inflatie'
+ ,'Uitvoer/inflatie'
+ ,'Inflatie')
+ ,1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('Invoer/inflatie','Uitvoer/inflatie','Inflatie'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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
Invoer/inflatie Uitvoer/inflatie Inflatie M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 10284.500 1.0383514 1.4 1 0 0 0 0 0 0 0 0 0 0
2 12792.000 0.9330311 1.3 0 1 0 0 0 0 0 0 0 0 0
3 12823.615 0.9327831 1.3 0 0 1 0 0 0 0 0 0 0 0
4 13845.667 0.9537554 1.2 0 0 0 1 0 0 0 0 0 0 0
5 15335.636 1.0098657 1.1 0 0 0 0 1 0 0 0 0 0 0
6 11188.500 0.9795325 1.4 0 0 0 0 0 1 0 0 0 0 0
7 13633.250 0.9865108 1.2 0 0 0 0 0 0 1 0 0 0 0
8 12298.467 0.9646613 1.5 0 0 0 0 0 0 0 1 0 0 0
9 15353.636 0.9467618 1.1 0 0 0 0 0 0 0 0 1 0 0
10 12696.154 0.9590689 1.3 0 0 0 0 0 0 0 0 0 1 0
11 12213.933 0.9857101 1.5 0 0 0 0 0 0 0 0 0 0 1
12 13683.727 0.9258216 1.1 0 0 0 0 0 0 0 0 0 0 0
13 11214.143 1.0368653 1.4 1 0 0 0 0 0 0 0 0 0 0
14 13950.231 0.9444436 1.3 0 1 0 0 0 0 0 0 0 0 0
15 11179.133 0.9449018 1.5 0 0 1 0 0 0 0 0 0 0 0
16 11801.875 0.9891514 1.6 0 0 0 1 0 0 0 0 0 0 0
17 11188.824 1.0543616 1.7 0 0 0 0 1 0 0 0 0 0 0
18 16456.273 1.0334789 1.1 0 0 0 0 0 1 0 0 0 0 0
19 11110.062 1.0013689 1.6 0 0 0 0 0 0 1 0 0 0 0
20 16530.692 1.0198126 1.3 0 0 0 0 0 0 0 1 0 0 0
21 10038.412 0.9939022 1.7 0 0 0 0 0 0 0 0 1 0 0
22 11681.250 0.9614445 1.6 0 0 0 0 0 0 0 0 0 1 0
23 11148.882 0.9574497 1.7 0 0 0 0 0 0 0 0 0 0 1
24 8631.000 0.9330864 1.9 0 0 0 0 0 0 0 0 0 0 0
25 9386.444 1.0451705 1.8 1 0 0 0 0 0 0 0 0 0 0
26 9764.737 0.9531663 1.9 0 1 0 0 0 0 0 0 0 0 0
27 12043.750 0.9667822 1.6 0 0 1 0 0 0 0 0 0 0 0
28 12948.067 0.9729926 1.5 0 0 0 1 0 0 0 0 0 0 0
29 10987.125 1.0136075 1.6 0 0 0 0 1 0 0 0 0 0 0
30 11648.312 0.9848395 1.6 0 0 0 0 0 1 0 0 0 0 0
31 10633.353 0.9732208 1.7 0 0 0 0 0 0 1 0 0 0 0
32 10219.300 0.9572846 2.0 0 0 0 0 0 0 0 1 0 0 0
33 9037.600 0.9720672 2.0 0 0 0 0 0 0 0 0 1 0 0
34 10296.316 0.9868789 1.9 0 0 0 0 0 0 0 0 0 1 0
35 11705.412 0.9546545 1.7 0 0 0 0 0 0 0 0 0 0 1
36 10681.944 0.9789870 1.8 0 0 0 0 0 0 0 0 0 0 0
37 9362.947 1.0030560 1.9 1 0 0 0 0 0 0 0 0 0 0
38 11306.353 0.9700812 1.7 0 1 0 0 0 0 0 0 0 0 0
39 10984.450 0.9913764 2.0 0 0 1 0 0 0 0 0 0 0 0
40 10062.619 1.0226090 2.1 0 0 0 1 0 0 0 0 0 0 0
41 8118.583 1.0899012 2.4 0 0 0 0 1 0 0 0 0 0 0
42 8867.480 1.0603736 2.5 0 0 0 0 0 1 0 0 0 0 0
43 8346.720 0.9859526 2.5 0 0 0 0 0 0 1 0 0 0 0
44 8529.308 1.0375122 2.6 0 0 0 0 0 0 0 1 0 0 0
45 10697.182 1.0253352 2.2 0 0 0 0 0 0 0 0 1 0 0
46 8591.840 1.0063766 2.5 0 0 0 0 0 0 0 0 0 1 0
47 8695.607 1.0187621 2.8 0 0 0 0 0 0 0 0 0 0 1
48 8125.571 1.0160185 2.8 0 0 0 0 0 0 0 0 0 0 0
49 7009.759 1.1124105 2.9 1 0 0 0 0 0 0 0 0 0 0
50 7883.467 1.0379037 3.0 0 1 0 0 0 0 0 0 0 0 0
51 7527.645 1.0454360 3.1 0 0 1 0 0 0 0 0 0 0 0
52 6763.759 1.0993543 2.9 0 0 0 1 0 0 0 0 0 0 0
53 6682.333 1.1019208 2.7 0 0 0 0 1 0 0 0 0 0 0
54 7855.682 1.0805750 2.2 0 0 0 0 0 1 0 0 0 0 0
55 6738.880 1.0243888 2.5 0 0 0 0 0 0 1 0 0 0 0
56 7895.435 1.0242822 2.3 0 0 0 0 0 0 0 1 0 0 0
57 6361.885 0.9938653 2.6 0 0 0 0 0 0 0 0 1 0 0
58 6935.957 0.9849352 2.3 0 0 0 0 0 0 0 0 0 1 0
59 8344.455 1.0057911 2.2 0 0 0 0 0 0 0 0 0 0 1
60 9107.944 0.9474283 1.8 0 0 0 0 0 0 0 0 0 0 0
t
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
19 19
20 20
21 21
22 22
23 23
24 24
25 25
26 26
27 27
28 28
29 29
30 30
31 31
32 32
33 33
34 34
35 35
36 36
37 37
38 38
39 39
40 40
41 41
42 42
43 43
44 44
45 45
46 46
47 47
48 48
49 49
50 50
51 51
52 52
53 53
54 54
55 55
56 56
57 57
58 58
59 59
60 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Uitvoer/inflatie` Inflatie M1
54.47 19736.99 -4260.49 -2600.41
M2 M3 M4 M5
511.42 397.47 -191.94 -1531.97
M6 M7 M8 M9
-844.71 -671.98 400.30 -172.66
M10 M11 t
-272.50 297.50 -26.43
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1881.0 -648.7 -116.6 659.2 2015.1
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 54.47 6051.68 0.009 0.99286
`Uitvoer/inflatie` 19736.99 6944.10 2.842 0.00671 **
Inflatie -4260.49 623.66 -6.831 1.81e-08 ***
M1 -2600.41 935.53 -2.780 0.00792 **
M2 511.42 707.68 0.723 0.47362
M3 397.47 714.55 0.556 0.58080
M4 -191.94 776.84 -0.247 0.80597
M5 -1531.97 952.05 -1.609 0.11458
M6 -844.71 859.82 -0.982 0.33114
M7 -671.98 732.05 -0.918 0.36354
M8 400.30 742.54 0.539 0.59248
M9 -172.66 712.22 -0.242 0.80956
M10 -272.50 701.30 -0.389 0.69943
M11 297.50 704.68 0.422 0.67490
t -26.43 16.56 -1.596 0.11743
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1089 on 45 degrees of freedom
Multiple R-squared: 0.8522, Adjusted R-squared: 0.8062
F-statistic: 18.54 on 14 and 45 DF, p-value: 3.516e-14
> 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,] 0.039655797 0.07931159 0.96034420
[2,] 0.014973465 0.02994693 0.98502653
[3,] 0.007917292 0.01583458 0.99208271
[4,] 0.027772437 0.05554487 0.97222756
[5,] 0.038155460 0.07631092 0.96184454
[6,] 0.019259276 0.03851855 0.98074072
[7,] 0.160672905 0.32134581 0.83932709
[8,] 0.160048733 0.32009747 0.83995127
[9,] 0.238079351 0.47615870 0.76192065
[10,] 0.192011247 0.38402249 0.80798875
[11,] 0.244146068 0.48829214 0.75585393
[12,] 0.430489682 0.86097936 0.56951032
[13,] 0.359962259 0.71992452 0.64003774
[14,] 0.268688113 0.53737623 0.73131189
[15,] 0.208310107 0.41662021 0.79168989
[16,] 0.448381637 0.89676327 0.55161836
[17,] 0.385900411 0.77180082 0.61409959
[18,] 0.290994610 0.58198922 0.70900539
[19,] 0.464500039 0.92900008 0.53549996
[20,] 0.465150217 0.93030043 0.53484978
[21,] 0.590776866 0.81844627 0.40922313
[22,] 0.690864591 0.61827082 0.30913541
[23,] 0.563964122 0.87207176 0.43603588
[24,] 0.942904312 0.11419138 0.05709569
[25,] 0.919998620 0.16000276 0.08000138
> postscript(file="/var/www/html/rcomp/tmp/1bcl61258728080.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/2e6dl1258728080.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/3e0xq1258728080.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/4wtnf1258728080.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/5yj6x1258728080.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 = 60
Frequency = 1
1 2 3 4 5 6
-1672.36971 -597.60896 -420.71863 377.19232 1700.12493 -1231.00506
7 8 9 10 11 12
77.61704 -593.62753 1710.01676 -211.99419 -911.50892 360.03978
13 14 15 16 17 18
-396.22440 652.54539 -1135.11762 -343.84377 -451.43848 2011.05187
19 20 21 22 23 24
-717.45713 2015.14949 -1662.15051 321.53328 -249.51626 -1110.50897
25 26 27 28 29 30
-366.47535 -831.64239 40.86599 1012.39705 42.34936 610.50360
31 32 33 34 35 36
104.61247 237.38782 -636.68666 29.91790 679.35381 -74.38157
37 38 39 40 41 42
1184.46171 -158.80212 517.51954 21.13608 -606.43710 490.46923
43 44 45 46 47 48
1292.25490 -162.58989 1140.81485 814.08182 1407.96978 1216.01713
49 50 51 52 53 54
1250.60775 935.50808 997.45072 -1066.88168 -684.59871 -1881.01964
55 56 57 58 59 60
-757.02728 -1496.31989 -551.99443 -953.53881 -926.29842 -391.16637
> postscript(file="/var/www/html/rcomp/tmp/63w431258728080.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -1672.36971 NA
1 -597.60896 -1672.36971
2 -420.71863 -597.60896
3 377.19232 -420.71863
4 1700.12493 377.19232
5 -1231.00506 1700.12493
6 77.61704 -1231.00506
7 -593.62753 77.61704
8 1710.01676 -593.62753
9 -211.99419 1710.01676
10 -911.50892 -211.99419
11 360.03978 -911.50892
12 -396.22440 360.03978
13 652.54539 -396.22440
14 -1135.11762 652.54539
15 -343.84377 -1135.11762
16 -451.43848 -343.84377
17 2011.05187 -451.43848
18 -717.45713 2011.05187
19 2015.14949 -717.45713
20 -1662.15051 2015.14949
21 321.53328 -1662.15051
22 -249.51626 321.53328
23 -1110.50897 -249.51626
24 -366.47535 -1110.50897
25 -831.64239 -366.47535
26 40.86599 -831.64239
27 1012.39705 40.86599
28 42.34936 1012.39705
29 610.50360 42.34936
30 104.61247 610.50360
31 237.38782 104.61247
32 -636.68666 237.38782
33 29.91790 -636.68666
34 679.35381 29.91790
35 -74.38157 679.35381
36 1184.46171 -74.38157
37 -158.80212 1184.46171
38 517.51954 -158.80212
39 21.13608 517.51954
40 -606.43710 21.13608
41 490.46923 -606.43710
42 1292.25490 490.46923
43 -162.58989 1292.25490
44 1140.81485 -162.58989
45 814.08182 1140.81485
46 1407.96978 814.08182
47 1216.01713 1407.96978
48 1250.60775 1216.01713
49 935.50808 1250.60775
50 997.45072 935.50808
51 -1066.88168 997.45072
52 -684.59871 -1066.88168
53 -1881.01964 -684.59871
54 -757.02728 -1881.01964
55 -1496.31989 -757.02728
56 -551.99443 -1496.31989
57 -953.53881 -551.99443
58 -926.29842 -953.53881
59 -391.16637 -926.29842
60 NA -391.16637
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -597.60896 -1672.36971
[2,] -420.71863 -597.60896
[3,] 377.19232 -420.71863
[4,] 1700.12493 377.19232
[5,] -1231.00506 1700.12493
[6,] 77.61704 -1231.00506
[7,] -593.62753 77.61704
[8,] 1710.01676 -593.62753
[9,] -211.99419 1710.01676
[10,] -911.50892 -211.99419
[11,] 360.03978 -911.50892
[12,] -396.22440 360.03978
[13,] 652.54539 -396.22440
[14,] -1135.11762 652.54539
[15,] -343.84377 -1135.11762
[16,] -451.43848 -343.84377
[17,] 2011.05187 -451.43848
[18,] -717.45713 2011.05187
[19,] 2015.14949 -717.45713
[20,] -1662.15051 2015.14949
[21,] 321.53328 -1662.15051
[22,] -249.51626 321.53328
[23,] -1110.50897 -249.51626
[24,] -366.47535 -1110.50897
[25,] -831.64239 -366.47535
[26,] 40.86599 -831.64239
[27,] 1012.39705 40.86599
[28,] 42.34936 1012.39705
[29,] 610.50360 42.34936
[30,] 104.61247 610.50360
[31,] 237.38782 104.61247
[32,] -636.68666 237.38782
[33,] 29.91790 -636.68666
[34,] 679.35381 29.91790
[35,] -74.38157 679.35381
[36,] 1184.46171 -74.38157
[37,] -158.80212 1184.46171
[38,] 517.51954 -158.80212
[39,] 21.13608 517.51954
[40,] -606.43710 21.13608
[41,] 490.46923 -606.43710
[42,] 1292.25490 490.46923
[43,] -162.58989 1292.25490
[44,] 1140.81485 -162.58989
[45,] 814.08182 1140.81485
[46,] 1407.96978 814.08182
[47,] 1216.01713 1407.96978
[48,] 1250.60775 1216.01713
[49,] 935.50808 1250.60775
[50,] 997.45072 935.50808
[51,] -1066.88168 997.45072
[52,] -684.59871 -1066.88168
[53,] -1881.01964 -684.59871
[54,] -757.02728 -1881.01964
[55,] -1496.31989 -757.02728
[56,] -551.99443 -1496.31989
[57,] -953.53881 -551.99443
[58,] -926.29842 -953.53881
[59,] -391.16637 -926.29842
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -597.60896 -1672.36971
2 -420.71863 -597.60896
3 377.19232 -420.71863
4 1700.12493 377.19232
5 -1231.00506 1700.12493
6 77.61704 -1231.00506
7 -593.62753 77.61704
8 1710.01676 -593.62753
9 -211.99419 1710.01676
10 -911.50892 -211.99419
11 360.03978 -911.50892
12 -396.22440 360.03978
13 652.54539 -396.22440
14 -1135.11762 652.54539
15 -343.84377 -1135.11762
16 -451.43848 -343.84377
17 2011.05187 -451.43848
18 -717.45713 2011.05187
19 2015.14949 -717.45713
20 -1662.15051 2015.14949
21 321.53328 -1662.15051
22 -249.51626 321.53328
23 -1110.50897 -249.51626
24 -366.47535 -1110.50897
25 -831.64239 -366.47535
26 40.86599 -831.64239
27 1012.39705 40.86599
28 42.34936 1012.39705
29 610.50360 42.34936
30 104.61247 610.50360
31 237.38782 104.61247
32 -636.68666 237.38782
33 29.91790 -636.68666
34 679.35381 29.91790
35 -74.38157 679.35381
36 1184.46171 -74.38157
37 -158.80212 1184.46171
38 517.51954 -158.80212
39 21.13608 517.51954
40 -606.43710 21.13608
41 490.46923 -606.43710
42 1292.25490 490.46923
43 -162.58989 1292.25490
44 1140.81485 -162.58989
45 814.08182 1140.81485
46 1407.96978 814.08182
47 1216.01713 1407.96978
48 1250.60775 1216.01713
49 935.50808 1250.60775
50 997.45072 935.50808
51 -1066.88168 997.45072
52 -684.59871 -1066.88168
53 -1881.01964 -684.59871
54 -757.02728 -1881.01964
55 -1496.31989 -757.02728
56 -551.99443 -1496.31989
57 -953.53881 -551.99443
58 -926.29842 -953.53881
59 -391.16637 -926.29842
> 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/729ga1258728080.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/83jlm1258728080.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/9v68j1258728080.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/10qj7y1258728080.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/11lp8a1258728080.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/12ywwt1258728080.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/13bx7d1258728080.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/148a1q1258728080.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/15woz81258728080.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/16trp61258728080.tab")
+ }
>
> system("convert tmp/1bcl61258728080.ps tmp/1bcl61258728080.png")
> system("convert tmp/2e6dl1258728080.ps tmp/2e6dl1258728080.png")
> system("convert tmp/3e0xq1258728080.ps tmp/3e0xq1258728080.png")
> system("convert tmp/4wtnf1258728080.ps tmp/4wtnf1258728080.png")
> system("convert tmp/5yj6x1258728080.ps tmp/5yj6x1258728080.png")
> system("convert tmp/63w431258728080.ps tmp/63w431258728080.png")
> system("convert tmp/729ga1258728080.ps tmp/729ga1258728080.png")
> system("convert tmp/83jlm1258728080.ps tmp/83jlm1258728080.png")
> system("convert tmp/9v68j1258728080.ps tmp/9v68j1258728080.png")
> system("convert tmp/10qj7y1258728080.ps tmp/10qj7y1258728080.png")
>
>
> proc.time()
user system elapsed
2.434 1.564 2.909