R version 2.7.0 (2008-04-22)
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(180144
+ ,8955.5
+ ,173666
+ ,10423.9
+ ,165688
+ ,11617.2
+ ,161570
+ ,9391.1
+ ,156145
+ ,10872
+ ,153730
+ ,10230.4
+ ,182698
+ ,9221
+ ,200765
+ ,9428.6
+ ,176512
+ ,10934.5
+ ,166618
+ ,10986
+ ,158644
+ ,11724.6
+ ,159585
+ ,11180.9
+ ,163095
+ ,11163.2
+ ,159044
+ ,11240.9
+ ,155511
+ ,12107.1
+ ,153745
+ ,10762.3
+ ,150569
+ ,11340.4
+ ,150605
+ ,11266.8
+ ,179612
+ ,9542.7
+ ,194690
+ ,9227.7
+ ,189917
+ ,10571.9
+ ,184128
+ ,10774.4
+ ,175335
+ ,10392.8
+ ,179566
+ ,9920.2
+ ,181140
+ ,9884.9
+ ,177876
+ ,10174.5
+ ,175041
+ ,11395.4
+ ,169292
+ ,10760.2
+ ,166070
+ ,10570.1
+ ,166972
+ ,10536
+ ,206348
+ ,9902.6
+ ,215706
+ ,8889
+ ,202108
+ ,10837.3
+ ,195411
+ ,11624.1
+ ,193111
+ ,10509
+ ,195198
+ ,10984.9
+ ,198770
+ ,10649.1
+ ,194163
+ ,10855.7
+ ,190420
+ ,11677.4
+ ,189733
+ ,10760.2
+ ,186029
+ ,10046.2
+ ,191531
+ ,10772.8
+ ,232571
+ ,9987.7
+ ,243477
+ ,8638.7
+ ,227247
+ ,11063.7
+ ,217859
+ ,11855.7
+ ,208679
+ ,10684.5
+ ,213188
+ ,11337.4
+ ,216234
+ ,10478
+ ,213586
+ ,11123.9
+ ,209465
+ ,12909.3
+ ,204045
+ ,11339.9
+ ,200237
+ ,10462.2
+ ,203666
+ ,12733.5
+ ,241476
+ ,10519.2
+ ,260307
+ ,10414.9
+ ,243324
+ ,12476.8
+ ,244460
+ ,12384.6
+ ,233575
+ ,12266.7
+ ,237217
+ ,12919.9
+ ,235243
+ ,11497.3
+ ,230354
+ ,12142
+ ,227184
+ ,13919.4
+ ,221678
+ ,12656.8
+ ,217142
+ ,12034.1
+ ,219452
+ ,13199.7
+ ,256446
+ ,10881.3
+ ,265845
+ ,11301.2
+ ,248624
+ ,13643.9
+ ,241114
+ ,12517
+ ,229245
+ ,13981.1
+ ,231805
+ ,14275.7
+ ,219277
+ ,13435
+ ,219313
+ ,13565.7
+ ,212610
+ ,16216.3
+ ,214771
+ ,12970
+ ,211142
+ ,14079.9
+ ,211457
+ ,14235
+ ,240048
+ ,12213.4
+ ,240636
+ ,12581
+ ,230580
+ ,14130.4
+ ,208795
+ ,14210.8
+ ,197922
+ ,14378.5
+ ,194596
+ ,13142.8)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'Europa')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','Europa'),1:84))
> 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
werkloosheid Europa M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 180144 8955.5 1 0 0 0 0 0 0 0 0 0 0 1
2 173666 10423.9 0 1 0 0 0 0 0 0 0 0 0 2
3 165688 11617.2 0 0 1 0 0 0 0 0 0 0 0 3
4 161570 9391.1 0 0 0 1 0 0 0 0 0 0 0 4
5 156145 10872.0 0 0 0 0 1 0 0 0 0 0 0 5
6 153730 10230.4 0 0 0 0 0 1 0 0 0 0 0 6
7 182698 9221.0 0 0 0 0 0 0 1 0 0 0 0 7
8 200765 9428.6 0 0 0 0 0 0 0 1 0 0 0 8
9 176512 10934.5 0 0 0 0 0 0 0 0 1 0 0 9
10 166618 10986.0 0 0 0 0 0 0 0 0 0 1 0 10
11 158644 11724.6 0 0 0 0 0 0 0 0 0 0 1 11
12 159585 11180.9 0 0 0 0 0 0 0 0 0 0 0 12
13 163095 11163.2 1 0 0 0 0 0 0 0 0 0 0 13
14 159044 11240.9 0 1 0 0 0 0 0 0 0 0 0 14
15 155511 12107.1 0 0 1 0 0 0 0 0 0 0 0 15
16 153745 10762.3 0 0 0 1 0 0 0 0 0 0 0 16
17 150569 11340.4 0 0 0 0 1 0 0 0 0 0 0 17
18 150605 11266.8 0 0 0 0 0 1 0 0 0 0 0 18
19 179612 9542.7 0 0 0 0 0 0 1 0 0 0 0 19
20 194690 9227.7 0 0 0 0 0 0 0 1 0 0 0 20
21 189917 10571.9 0 0 0 0 0 0 0 0 1 0 0 21
22 184128 10774.4 0 0 0 0 0 0 0 0 0 1 0 22
23 175335 10392.8 0 0 0 0 0 0 0 0 0 0 1 23
24 179566 9920.2 0 0 0 0 0 0 0 0 0 0 0 24
25 181140 9884.9 1 0 0 0 0 0 0 0 0 0 0 25
26 177876 10174.5 0 1 0 0 0 0 0 0 0 0 0 26
27 175041 11395.4 0 0 1 0 0 0 0 0 0 0 0 27
28 169292 10760.2 0 0 0 1 0 0 0 0 0 0 0 28
29 166070 10570.1 0 0 0 0 1 0 0 0 0 0 0 29
30 166972 10536.0 0 0 0 0 0 1 0 0 0 0 0 30
31 206348 9902.6 0 0 0 0 0 0 1 0 0 0 0 31
32 215706 8889.0 0 0 0 0 0 0 0 1 0 0 0 32
33 202108 10837.3 0 0 0 0 0 0 0 0 1 0 0 33
34 195411 11624.1 0 0 0 0 0 0 0 0 0 1 0 34
35 193111 10509.0 0 0 0 0 0 0 0 0 0 0 1 35
36 195198 10984.9 0 0 0 0 0 0 0 0 0 0 0 36
37 198770 10649.1 1 0 0 0 0 0 0 0 0 0 0 37
38 194163 10855.7 0 1 0 0 0 0 0 0 0 0 0 38
39 190420 11677.4 0 0 1 0 0 0 0 0 0 0 0 39
40 189733 10760.2 0 0 0 1 0 0 0 0 0 0 0 40
41 186029 10046.2 0 0 0 0 1 0 0 0 0 0 0 41
42 191531 10772.8 0 0 0 0 0 1 0 0 0 0 0 42
43 232571 9987.7 0 0 0 0 0 0 1 0 0 0 0 43
44 243477 8638.7 0 0 0 0 0 0 0 1 0 0 0 44
45 227247 11063.7 0 0 0 0 0 0 0 0 1 0 0 45
46 217859 11855.7 0 0 0 0 0 0 0 0 0 1 0 46
47 208679 10684.5 0 0 0 0 0 0 0 0 0 0 1 47
48 213188 11337.4 0 0 0 0 0 0 0 0 0 0 0 48
49 216234 10478.0 1 0 0 0 0 0 0 0 0 0 0 49
50 213586 11123.9 0 1 0 0 0 0 0 0 0 0 0 50
51 209465 12909.3 0 0 1 0 0 0 0 0 0 0 0 51
52 204045 11339.9 0 0 0 1 0 0 0 0 0 0 0 52
53 200237 10462.2 0 0 0 0 1 0 0 0 0 0 0 53
54 203666 12733.5 0 0 0 0 0 1 0 0 0 0 0 54
55 241476 10519.2 0 0 0 0 0 0 1 0 0 0 0 55
56 260307 10414.9 0 0 0 0 0 0 0 1 0 0 0 56
57 243324 12476.8 0 0 0 0 0 0 0 0 1 0 0 57
58 244460 12384.6 0 0 0 0 0 0 0 0 0 1 0 58
59 233575 12266.7 0 0 0 0 0 0 0 0 0 0 1 59
60 237217 12919.9 0 0 0 0 0 0 0 0 0 0 0 60
61 235243 11497.3 1 0 0 0 0 0 0 0 0 0 0 61
62 230354 12142.0 0 1 0 0 0 0 0 0 0 0 0 62
63 227184 13919.4 0 0 1 0 0 0 0 0 0 0 0 63
64 221678 12656.8 0 0 0 1 0 0 0 0 0 0 0 64
65 217142 12034.1 0 0 0 0 1 0 0 0 0 0 0 65
66 219452 13199.7 0 0 0 0 0 1 0 0 0 0 0 66
67 256446 10881.3 0 0 0 0 0 0 1 0 0 0 0 67
68 265845 11301.2 0 0 0 0 0 0 0 1 0 0 0 68
69 248624 13643.9 0 0 0 0 0 0 0 0 1 0 0 69
70 241114 12517.0 0 0 0 0 0 0 0 0 0 1 0 70
71 229245 13981.1 0 0 0 0 0 0 0 0 0 0 1 71
72 231805 14275.7 0 0 0 0 0 0 0 0 0 0 0 72
73 219277 13435.0 1 0 0 0 0 0 0 0 0 0 0 73
74 219313 13565.7 0 1 0 0 0 0 0 0 0 0 0 74
75 212610 16216.3 0 0 1 0 0 0 0 0 0 0 0 75
76 214771 12970.0 0 0 0 1 0 0 0 0 0 0 0 76
77 211142 14079.9 0 0 0 0 1 0 0 0 0 0 0 77
78 211457 14235.0 0 0 0 0 0 1 0 0 0 0 0 78
79 240048 12213.4 0 0 0 0 0 0 1 0 0 0 0 79
80 240636 12581.0 0 0 0 0 0 0 0 1 0 0 0 80
81 230580 14130.4 0 0 0 0 0 0 0 0 1 0 0 81
82 208795 14210.8 0 0 0 0 0 0 0 0 0 1 0 82
83 197922 14378.5 0 0 0 0 0 0 0 0 0 0 1 83
84 194596 13142.8 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Europa M1 M2 M3 M4
217808.739 -6.194 3991.801 2150.081 5487.865 -8641.811
M5 M6 M7 M8 M9 M10
-13100.247 -9708.443 14153.128 23112.644 18835.388 9682.321
M11 t
-730.588 1206.215
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-43132.4 -7057.3 -875.6 8768.1 27057.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 217808.739 20018.016 10.881 < 2e-16 ***
Europa -6.194 1.960 -3.160 0.00233 **
M1 3991.801 7028.935 0.568 0.57191
M2 2150.081 6930.068 0.310 0.75729
M3 5487.865 7356.344 0.746 0.45816
M4 -8641.811 6952.601 -1.243 0.21803
M5 -13100.247 6936.887 -1.888 0.06310 .
M6 -9708.443 6914.565 -1.404 0.16472
M7 14153.128 7442.106 1.902 0.06132 .
M8 23112.644 7679.229 3.010 0.00363 **
M9 18835.388 6904.916 2.728 0.00805 **
M10 9682.321 6908.108 1.402 0.16546
M11 -730.588 6900.239 -0.106 0.91598
t 1206.215 106.179 11.360 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12910 on 70 degrees of freedom
Multiple R-squared: 0.838, Adjusted R-squared: 0.808
F-statistic: 27.86 on 13 and 70 DF, p-value: < 2.2e-16
> 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.338671e-02 2.677342e-02 0.9866133
[2,] 8.273688e-03 1.654738e-02 0.9917263
[3,] 2.283840e-03 4.567680e-03 0.9977162
[4,] 5.713794e-04 1.142759e-03 0.9994286
[5,] 2.671137e-03 5.342274e-03 0.9973289
[6,] 5.212029e-03 1.042406e-02 0.9947880
[7,] 1.911605e-03 3.823210e-03 0.9980884
[8,] 7.549318e-04 1.509864e-03 0.9992451
[9,] 2.862731e-04 5.725462e-04 0.9997137
[10,] 9.917678e-05 1.983536e-04 0.9999008
[11,] 4.162076e-05 8.324153e-05 0.9999584
[12,] 9.260572e-05 1.852114e-04 0.9999074
[13,] 3.559712e-05 7.119425e-05 0.9999644
[14,] 2.111543e-05 4.223087e-05 0.9999789
[15,] 4.352951e-04 8.705903e-04 0.9995647
[16,] 2.738751e-04 5.477501e-04 0.9997261
[17,] 2.781644e-04 5.563287e-04 0.9997218
[18,] 6.740742e-04 1.348148e-03 0.9993259
[19,] 5.386490e-04 1.077298e-03 0.9994614
[20,] 1.063909e-03 2.127817e-03 0.9989361
[21,] 1.282679e-03 2.565358e-03 0.9987173
[22,] 1.314140e-03 2.628280e-03 0.9986859
[23,] 1.023438e-03 2.046875e-03 0.9989766
[24,] 1.545763e-03 3.091526e-03 0.9984542
[25,] 1.350520e-03 2.701039e-03 0.9986495
[26,] 2.020498e-03 4.040995e-03 0.9979795
[27,] 9.503170e-03 1.900634e-02 0.9904968
[28,] 8.711735e-03 1.742347e-02 0.9912883
[29,] 1.151807e-02 2.303613e-02 0.9884819
[30,] 2.315877e-02 4.631755e-02 0.9768412
[31,] 1.869741e-02 3.739483e-02 0.9813026
[32,] 2.212698e-02 4.425396e-02 0.9778730
[33,] 1.913822e-02 3.827644e-02 0.9808618
[34,] 2.019622e-02 4.039245e-02 0.9798038
[35,] 2.316710e-02 4.633420e-02 0.9768329
[36,] 4.170162e-02 8.340324e-02 0.9582984
[37,] 6.932659e-02 1.386532e-01 0.9306734
[38,] 2.235918e-01 4.471836e-01 0.7764082
[39,] 5.281850e-01 9.436300e-01 0.4718150
[40,] 6.255185e-01 7.489629e-01 0.3744815
[41,] 8.367499e-01 3.265003e-01 0.1632501
[42,] 8.571205e-01 2.857590e-01 0.1428795
[43,] 8.298728e-01 3.402544e-01 0.1701272
[44,] 8.151056e-01 3.697888e-01 0.1848944
[45,] 7.365479e-01 5.269043e-01 0.2634521
[46,] 6.603495e-01 6.793010e-01 0.3396505
[47,] 5.509328e-01 8.981344e-01 0.4490672
[48,] 6.698013e-01 6.603975e-01 0.3301987
[49,] 6.232879e-01 7.534242e-01 0.3767121
[50,] 7.649516e-01 4.700968e-01 0.2350484
[51,] 7.287244e-01 5.425511e-01 0.2712756
> postscript(file="/var/www/html/rcomp/tmp/1tem81229725606.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/2qegv1229725606.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/3ehz91229725606.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/4c8hs1229725606.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/5su4q1229725606.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 = 84
Frequency = 1
1 2 3 4 5 6
12604.8256 15857.1440 10726.0733 5743.7632 12743.2178 1756.3278
7 8 9 10 11 12
-595.3687 8591.7102 -3263.1731 -4891.3467 916.0061 -3447.3065
13 14 15 16 17 18
-5244.9516 -8179.1965 -10891.2229 -8063.0349 -4406.2427 -9424.1174
19 20 21 22 23 24
-16163.4433 -13202.1863 -6578.5892 -3166.5157 -5116.3333 -5749.2746
25 26 27 28 29 30
-9591.9287 -10426.7308 -10243.8557 -7003.6260 -8150.8257 -12058.0497
31 32 33 34 35 36
-1672.9191 -8758.5732 -7218.3685 -1095.3225 -1095.2111 2002.5617
37 38 39 40 41 42
-1703.2960 -4395.1743 -7592.8199 -1037.2104 -5911.2838 -506.9685
43 44 45 46 47 48
10602.5794 2987.5621 4848.2982 8312.5514 1085.1969 7701.2526
49 50 51 52 53 54
226.3794 2214.3887 4607.6012 2390.6867 -3601.2936 9297.4131
55 56 57 58 59 60
8324.9407 16344.2080 15203.0158 23714.8091 21306.2673 27057.1811
61 62 63 64 65 66
11074.0222 10813.5991 14108.2621 13705.5713 8564.9817 13496.3265
67 68 69 70 71 72
11063.0911 12897.0902 13257.0858 6714.2691 13120.1433 15568.0002
73 74 75 76 77 78
-7365.0510 -5884.0303 -714.0381 -5736.1497 761.4463 -2560.9317
79 80 81 82 83 84
-11558.8802 -18859.8110 -16248.2689 -29588.4448 -30216.0691 -43132.4145
> postscript(file="/var/www/html/rcomp/tmp/6jmgn1229725606.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 12604.8256 NA
1 15857.1440 12604.8256
2 10726.0733 15857.1440
3 5743.7632 10726.0733
4 12743.2178 5743.7632
5 1756.3278 12743.2178
6 -595.3687 1756.3278
7 8591.7102 -595.3687
8 -3263.1731 8591.7102
9 -4891.3467 -3263.1731
10 916.0061 -4891.3467
11 -3447.3065 916.0061
12 -5244.9516 -3447.3065
13 -8179.1965 -5244.9516
14 -10891.2229 -8179.1965
15 -8063.0349 -10891.2229
16 -4406.2427 -8063.0349
17 -9424.1174 -4406.2427
18 -16163.4433 -9424.1174
19 -13202.1863 -16163.4433
20 -6578.5892 -13202.1863
21 -3166.5157 -6578.5892
22 -5116.3333 -3166.5157
23 -5749.2746 -5116.3333
24 -9591.9287 -5749.2746
25 -10426.7308 -9591.9287
26 -10243.8557 -10426.7308
27 -7003.6260 -10243.8557
28 -8150.8257 -7003.6260
29 -12058.0497 -8150.8257
30 -1672.9191 -12058.0497
31 -8758.5732 -1672.9191
32 -7218.3685 -8758.5732
33 -1095.3225 -7218.3685
34 -1095.2111 -1095.3225
35 2002.5617 -1095.2111
36 -1703.2960 2002.5617
37 -4395.1743 -1703.2960
38 -7592.8199 -4395.1743
39 -1037.2104 -7592.8199
40 -5911.2838 -1037.2104
41 -506.9685 -5911.2838
42 10602.5794 -506.9685
43 2987.5621 10602.5794
44 4848.2982 2987.5621
45 8312.5514 4848.2982
46 1085.1969 8312.5514
47 7701.2526 1085.1969
48 226.3794 7701.2526
49 2214.3887 226.3794
50 4607.6012 2214.3887
51 2390.6867 4607.6012
52 -3601.2936 2390.6867
53 9297.4131 -3601.2936
54 8324.9407 9297.4131
55 16344.2080 8324.9407
56 15203.0158 16344.2080
57 23714.8091 15203.0158
58 21306.2673 23714.8091
59 27057.1811 21306.2673
60 11074.0222 27057.1811
61 10813.5991 11074.0222
62 14108.2621 10813.5991
63 13705.5713 14108.2621
64 8564.9817 13705.5713
65 13496.3265 8564.9817
66 11063.0911 13496.3265
67 12897.0902 11063.0911
68 13257.0858 12897.0902
69 6714.2691 13257.0858
70 13120.1433 6714.2691
71 15568.0002 13120.1433
72 -7365.0510 15568.0002
73 -5884.0303 -7365.0510
74 -714.0381 -5884.0303
75 -5736.1497 -714.0381
76 761.4463 -5736.1497
77 -2560.9317 761.4463
78 -11558.8802 -2560.9317
79 -18859.8110 -11558.8802
80 -16248.2689 -18859.8110
81 -29588.4448 -16248.2689
82 -30216.0691 -29588.4448
83 -43132.4145 -30216.0691
84 NA -43132.4145
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 15857.1440 12604.8256
[2,] 10726.0733 15857.1440
[3,] 5743.7632 10726.0733
[4,] 12743.2178 5743.7632
[5,] 1756.3278 12743.2178
[6,] -595.3687 1756.3278
[7,] 8591.7102 -595.3687
[8,] -3263.1731 8591.7102
[9,] -4891.3467 -3263.1731
[10,] 916.0061 -4891.3467
[11,] -3447.3065 916.0061
[12,] -5244.9516 -3447.3065
[13,] -8179.1965 -5244.9516
[14,] -10891.2229 -8179.1965
[15,] -8063.0349 -10891.2229
[16,] -4406.2427 -8063.0349
[17,] -9424.1174 -4406.2427
[18,] -16163.4433 -9424.1174
[19,] -13202.1863 -16163.4433
[20,] -6578.5892 -13202.1863
[21,] -3166.5157 -6578.5892
[22,] -5116.3333 -3166.5157
[23,] -5749.2746 -5116.3333
[24,] -9591.9287 -5749.2746
[25,] -10426.7308 -9591.9287
[26,] -10243.8557 -10426.7308
[27,] -7003.6260 -10243.8557
[28,] -8150.8257 -7003.6260
[29,] -12058.0497 -8150.8257
[30,] -1672.9191 -12058.0497
[31,] -8758.5732 -1672.9191
[32,] -7218.3685 -8758.5732
[33,] -1095.3225 -7218.3685
[34,] -1095.2111 -1095.3225
[35,] 2002.5617 -1095.2111
[36,] -1703.2960 2002.5617
[37,] -4395.1743 -1703.2960
[38,] -7592.8199 -4395.1743
[39,] -1037.2104 -7592.8199
[40,] -5911.2838 -1037.2104
[41,] -506.9685 -5911.2838
[42,] 10602.5794 -506.9685
[43,] 2987.5621 10602.5794
[44,] 4848.2982 2987.5621
[45,] 8312.5514 4848.2982
[46,] 1085.1969 8312.5514
[47,] 7701.2526 1085.1969
[48,] 226.3794 7701.2526
[49,] 2214.3887 226.3794
[50,] 4607.6012 2214.3887
[51,] 2390.6867 4607.6012
[52,] -3601.2936 2390.6867
[53,] 9297.4131 -3601.2936
[54,] 8324.9407 9297.4131
[55,] 16344.2080 8324.9407
[56,] 15203.0158 16344.2080
[57,] 23714.8091 15203.0158
[58,] 21306.2673 23714.8091
[59,] 27057.1811 21306.2673
[60,] 11074.0222 27057.1811
[61,] 10813.5991 11074.0222
[62,] 14108.2621 10813.5991
[63,] 13705.5713 14108.2621
[64,] 8564.9817 13705.5713
[65,] 13496.3265 8564.9817
[66,] 11063.0911 13496.3265
[67,] 12897.0902 11063.0911
[68,] 13257.0858 12897.0902
[69,] 6714.2691 13257.0858
[70,] 13120.1433 6714.2691
[71,] 15568.0002 13120.1433
[72,] -7365.0510 15568.0002
[73,] -5884.0303 -7365.0510
[74,] -714.0381 -5884.0303
[75,] -5736.1497 -714.0381
[76,] 761.4463 -5736.1497
[77,] -2560.9317 761.4463
[78,] -11558.8802 -2560.9317
[79,] -18859.8110 -11558.8802
[80,] -16248.2689 -18859.8110
[81,] -29588.4448 -16248.2689
[82,] -30216.0691 -29588.4448
[83,] -43132.4145 -30216.0691
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 15857.1440 12604.8256
2 10726.0733 15857.1440
3 5743.7632 10726.0733
4 12743.2178 5743.7632
5 1756.3278 12743.2178
6 -595.3687 1756.3278
7 8591.7102 -595.3687
8 -3263.1731 8591.7102
9 -4891.3467 -3263.1731
10 916.0061 -4891.3467
11 -3447.3065 916.0061
12 -5244.9516 -3447.3065
13 -8179.1965 -5244.9516
14 -10891.2229 -8179.1965
15 -8063.0349 -10891.2229
16 -4406.2427 -8063.0349
17 -9424.1174 -4406.2427
18 -16163.4433 -9424.1174
19 -13202.1863 -16163.4433
20 -6578.5892 -13202.1863
21 -3166.5157 -6578.5892
22 -5116.3333 -3166.5157
23 -5749.2746 -5116.3333
24 -9591.9287 -5749.2746
25 -10426.7308 -9591.9287
26 -10243.8557 -10426.7308
27 -7003.6260 -10243.8557
28 -8150.8257 -7003.6260
29 -12058.0497 -8150.8257
30 -1672.9191 -12058.0497
31 -8758.5732 -1672.9191
32 -7218.3685 -8758.5732
33 -1095.3225 -7218.3685
34 -1095.2111 -1095.3225
35 2002.5617 -1095.2111
36 -1703.2960 2002.5617
37 -4395.1743 -1703.2960
38 -7592.8199 -4395.1743
39 -1037.2104 -7592.8199
40 -5911.2838 -1037.2104
41 -506.9685 -5911.2838
42 10602.5794 -506.9685
43 2987.5621 10602.5794
44 4848.2982 2987.5621
45 8312.5514 4848.2982
46 1085.1969 8312.5514
47 7701.2526 1085.1969
48 226.3794 7701.2526
49 2214.3887 226.3794
50 4607.6012 2214.3887
51 2390.6867 4607.6012
52 -3601.2936 2390.6867
53 9297.4131 -3601.2936
54 8324.9407 9297.4131
55 16344.2080 8324.9407
56 15203.0158 16344.2080
57 23714.8091 15203.0158
58 21306.2673 23714.8091
59 27057.1811 21306.2673
60 11074.0222 27057.1811
61 10813.5991 11074.0222
62 14108.2621 10813.5991
63 13705.5713 14108.2621
64 8564.9817 13705.5713
65 13496.3265 8564.9817
66 11063.0911 13496.3265
67 12897.0902 11063.0911
68 13257.0858 12897.0902
69 6714.2691 13257.0858
70 13120.1433 6714.2691
71 15568.0002 13120.1433
72 -7365.0510 15568.0002
73 -5884.0303 -7365.0510
74 -714.0381 -5884.0303
75 -5736.1497 -714.0381
76 761.4463 -5736.1497
77 -2560.9317 761.4463
78 -11558.8802 -2560.9317
79 -18859.8110 -11558.8802
80 -16248.2689 -18859.8110
81 -29588.4448 -16248.2689
82 -30216.0691 -29588.4448
83 -43132.4145 -30216.0691
> 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/79dhm1229725606.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/8v7sp1229725606.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/98n1p1229725606.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/10j9i21229725606.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/11ed8z1229725606.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/12z0un1229725606.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/13flyh1229725606.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/14dad91229725606.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/15djnu1229725606.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/16pq401229725606.tab")
+ }
>
> system("convert tmp/1tem81229725606.ps tmp/1tem81229725606.png")
> system("convert tmp/2qegv1229725606.ps tmp/2qegv1229725606.png")
> system("convert tmp/3ehz91229725606.ps tmp/3ehz91229725606.png")
> system("convert tmp/4c8hs1229725606.ps tmp/4c8hs1229725606.png")
> system("convert tmp/5su4q1229725606.ps tmp/5su4q1229725606.png")
> system("convert tmp/6jmgn1229725606.ps tmp/6jmgn1229725606.png")
> system("convert tmp/79dhm1229725606.ps tmp/79dhm1229725606.png")
> system("convert tmp/8v7sp1229725606.ps tmp/8v7sp1229725606.png")
> system("convert tmp/98n1p1229725606.ps tmp/98n1p1229725606.png")
> system("convert tmp/10j9i21229725606.ps tmp/10j9i21229725606.png")
>
>
> proc.time()
user system elapsed
5.586 2.814 5.939