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
+ ,1235.8
+ ,173666
+ ,1147.1
+ ,165688
+ ,1376.9
+ ,161570
+ ,1157.7
+ ,156145
+ ,1506
+ ,153730
+ ,1271.3
+ ,182698
+ ,1240.2
+ ,200765
+ ,1408.3
+ ,176512
+ ,1334.6
+ ,166618
+ ,1601.2
+ ,158644
+ ,1566.4
+ ,159585
+ ,1297.5
+ ,163095
+ ,1487.6
+ ,159044
+ ,1320.9
+ ,155511
+ ,1514
+ ,153745
+ ,1290.9
+ ,150569
+ ,1392.5
+ ,150605
+ ,1288.2
+ ,179612
+ ,1304.4
+ ,194690
+ ,1297.8
+ ,189917
+ ,1211
+ ,184128
+ ,1454
+ ,175335
+ ,1405.7
+ ,179566
+ ,1160.8
+ ,181140
+ ,1492.1
+ ,177876
+ ,1263
+ ,175041
+ ,1376.3
+ ,169292
+ ,1368.6
+ ,166070
+ ,1427.6
+ ,166972
+ ,1339.8
+ ,206348
+ ,1248.3
+ ,215706
+ ,1309.8
+ ,202108
+ ,1424
+ ,195411
+ ,1590.5
+ ,193111
+ ,1423.1
+ ,195198
+ ,1355.3
+ ,198770
+ ,1515
+ ,194163
+ ,1385.6
+ ,190420
+ ,1430
+ ,189733
+ ,1494.2
+ ,186029
+ ,1580.9
+ ,191531
+ ,1369.8
+ ,232571
+ ,1407.5
+ ,243477
+ ,1388.3
+ ,227247
+ ,1478.5
+ ,217859
+ ,1630.4
+ ,208679
+ ,1413.5
+ ,213188
+ ,1493.8
+ ,216234
+ ,1641.3
+ ,213586
+ ,1465
+ ,209465
+ ,1725.1
+ ,204045
+ ,1628.4
+ ,200237
+ ,1679.8
+ ,203666
+ ,1876
+ ,241476
+ ,1669.4
+ ,260307
+ ,1712.4
+ ,243324
+ ,1768.8
+ ,244460
+ ,1820.5
+ ,233575
+ ,1776.2
+ ,237217
+ ,1693.7
+ ,235243
+ ,1799.1
+ ,230354
+ ,1917.5
+ ,227184
+ ,1887.2
+ ,221678
+ ,1787.8
+ ,217142
+ ,1803.8
+ ,219452
+ ,2196.4
+ ,256446
+ ,1759.5
+ ,265845
+ ,2002.6
+ ,248624
+ ,2056.8
+ ,241114
+ ,1851.1
+ ,229245
+ ,1984.3
+ ,231805
+ ,1725.3
+ ,219277
+ ,2096.6
+ ,219313
+ ,1792.2
+ ,212610
+ ,2029.9
+ ,214771
+ ,1785.3
+ ,211142
+ ,2026.5
+ ,211457
+ ,1930.8
+ ,240048
+ ,1845.5
+ ,240636
+ ,1943.1
+ ,230580
+ ,2066.8
+ ,208795
+ ,2354.4
+ ,197922
+ ,2190.7
+ ,194596
+ ,1929.6)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'Azië')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','Azië'),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 Azi\353 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 180144 1235.8 1 0 0 0 0 0 0 0 0 0 0 1
2 173666 1147.1 0 1 0 0 0 0 0 0 0 0 0 2
3 165688 1376.9 0 0 1 0 0 0 0 0 0 0 0 3
4 161570 1157.7 0 0 0 1 0 0 0 0 0 0 0 4
5 156145 1506.0 0 0 0 0 1 0 0 0 0 0 0 5
6 153730 1271.3 0 0 0 0 0 1 0 0 0 0 0 6
7 182698 1240.2 0 0 0 0 0 0 1 0 0 0 0 7
8 200765 1408.3 0 0 0 0 0 0 0 1 0 0 0 8
9 176512 1334.6 0 0 0 0 0 0 0 0 1 0 0 9
10 166618 1601.2 0 0 0 0 0 0 0 0 0 1 0 10
11 158644 1566.4 0 0 0 0 0 0 0 0 0 0 1 11
12 159585 1297.5 0 0 0 0 0 0 0 0 0 0 0 12
13 163095 1487.6 1 0 0 0 0 0 0 0 0 0 0 13
14 159044 1320.9 0 1 0 0 0 0 0 0 0 0 0 14
15 155511 1514.0 0 0 1 0 0 0 0 0 0 0 0 15
16 153745 1290.9 0 0 0 1 0 0 0 0 0 0 0 16
17 150569 1392.5 0 0 0 0 1 0 0 0 0 0 0 17
18 150605 1288.2 0 0 0 0 0 1 0 0 0 0 0 18
19 179612 1304.4 0 0 0 0 0 0 1 0 0 0 0 19
20 194690 1297.8 0 0 0 0 0 0 0 1 0 0 0 20
21 189917 1211.0 0 0 0 0 0 0 0 0 1 0 0 21
22 184128 1454.0 0 0 0 0 0 0 0 0 0 1 0 22
23 175335 1405.7 0 0 0 0 0 0 0 0 0 0 1 23
24 179566 1160.8 0 0 0 0 0 0 0 0 0 0 0 24
25 181140 1492.1 1 0 0 0 0 0 0 0 0 0 0 25
26 177876 1263.0 0 1 0 0 0 0 0 0 0 0 0 26
27 175041 1376.3 0 0 1 0 0 0 0 0 0 0 0 27
28 169292 1368.6 0 0 0 1 0 0 0 0 0 0 0 28
29 166070 1427.6 0 0 0 0 1 0 0 0 0 0 0 29
30 166972 1339.8 0 0 0 0 0 1 0 0 0 0 0 30
31 206348 1248.3 0 0 0 0 0 0 1 0 0 0 0 31
32 215706 1309.8 0 0 0 0 0 0 0 1 0 0 0 32
33 202108 1424.0 0 0 0 0 0 0 0 0 1 0 0 33
34 195411 1590.5 0 0 0 0 0 0 0 0 0 1 0 34
35 193111 1423.1 0 0 0 0 0 0 0 0 0 0 1 35
36 195198 1355.3 0 0 0 0 0 0 0 0 0 0 0 36
37 198770 1515.0 1 0 0 0 0 0 0 0 0 0 0 37
38 194163 1385.6 0 1 0 0 0 0 0 0 0 0 0 38
39 190420 1430.0 0 0 1 0 0 0 0 0 0 0 0 39
40 189733 1494.2 0 0 0 1 0 0 0 0 0 0 0 40
41 186029 1580.9 0 0 0 0 1 0 0 0 0 0 0 41
42 191531 1369.8 0 0 0 0 0 1 0 0 0 0 0 42
43 232571 1407.5 0 0 0 0 0 0 1 0 0 0 0 43
44 243477 1388.3 0 0 0 0 0 0 0 1 0 0 0 44
45 227247 1478.5 0 0 0 0 0 0 0 0 1 0 0 45
46 217859 1630.4 0 0 0 0 0 0 0 0 0 1 0 46
47 208679 1413.5 0 0 0 0 0 0 0 0 0 0 1 47
48 213188 1493.8 0 0 0 0 0 0 0 0 0 0 0 48
49 216234 1641.3 1 0 0 0 0 0 0 0 0 0 0 49
50 213586 1465.0 0 1 0 0 0 0 0 0 0 0 0 50
51 209465 1725.1 0 0 1 0 0 0 0 0 0 0 0 51
52 204045 1628.4 0 0 0 1 0 0 0 0 0 0 0 52
53 200237 1679.8 0 0 0 0 1 0 0 0 0 0 0 53
54 203666 1876.0 0 0 0 0 0 1 0 0 0 0 0 54
55 241476 1669.4 0 0 0 0 0 0 1 0 0 0 0 55
56 260307 1712.4 0 0 0 0 0 0 0 1 0 0 0 56
57 243324 1768.8 0 0 0 0 0 0 0 0 1 0 0 57
58 244460 1820.5 0 0 0 0 0 0 0 0 0 1 0 58
59 233575 1776.2 0 0 0 0 0 0 0 0 0 0 1 59
60 237217 1693.7 0 0 0 0 0 0 0 0 0 0 0 60
61 235243 1799.1 1 0 0 0 0 0 0 0 0 0 0 61
62 230354 1917.5 0 1 0 0 0 0 0 0 0 0 0 62
63 227184 1887.2 0 0 1 0 0 0 0 0 0 0 0 63
64 221678 1787.8 0 0 0 1 0 0 0 0 0 0 0 64
65 217142 1803.8 0 0 0 0 1 0 0 0 0 0 0 65
66 219452 2196.4 0 0 0 0 0 1 0 0 0 0 0 66
67 256446 1759.5 0 0 0 0 0 0 1 0 0 0 0 67
68 265845 2002.6 0 0 0 0 0 0 0 1 0 0 0 68
69 248624 2056.8 0 0 0 0 0 0 0 0 1 0 0 69
70 241114 1851.1 0 0 0 0 0 0 0 0 0 1 0 70
71 229245 1984.3 0 0 0 0 0 0 0 0 0 0 1 71
72 231805 1725.3 0 0 0 0 0 0 0 0 0 0 0 72
73 219277 2096.6 1 0 0 0 0 0 0 0 0 0 0 73
74 219313 1792.2 0 1 0 0 0 0 0 0 0 0 0 74
75 212610 2029.9 0 0 1 0 0 0 0 0 0 0 0 75
76 214771 1785.3 0 0 0 1 0 0 0 0 0 0 0 76
77 211142 2026.5 0 0 0 0 1 0 0 0 0 0 0 77
78 211457 1930.8 0 0 0 0 0 1 0 0 0 0 0 78
79 240048 1845.5 0 0 0 0 0 0 1 0 0 0 0 79
80 240636 1943.1 0 0 0 0 0 0 0 1 0 0 0 80
81 230580 2066.8 0 0 0 0 0 0 0 0 1 0 0 81
82 208795 2354.4 0 0 0 0 0 0 0 0 0 1 0 82
83 197922 2190.7 0 0 0 0 0 0 0 0 0 0 1 83
84 194596 1929.6 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) `Azi\353` M1 M2 M3 M4
176169.70 -18.35 11367.99 3996.51 1049.59 -5241.38
M5 M6 M7 M8 M9 M10
-7910.89 -7962.45 23375.63 35551.02 20438.14 13286.70
M11 t
1914.16 1111.75
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-39544.6 -7915.7 -331.3 9759.0 25428.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 176169.70 14280.03 12.337 < 2e-16 ***
`Azi\353` -18.35 12.51 -1.467 0.14692
M1 11367.99 7705.62 1.475 0.14462
M2 3996.51 7315.09 0.546 0.58657
M3 1049.59 7656.65 0.137 0.89136
M4 -5241.38 7319.10 -0.716 0.47630
M5 -7910.89 7614.43 -1.039 0.30241
M6 -7962.45 7506.66 -1.061 0.29246
M7 23375.63 7276.42 3.213 0.00199 **
M8 35551.02 7371.12 4.823 7.97e-06 ***
M9 20438.14 7440.31 2.747 0.00764 **
M10 13286.70 7936.01 1.674 0.09855 .
M11 1914.16 7560.99 0.253 0.80088
t 1111.75 140.53 7.911 2.66e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13590 on 70 degrees of freedom
Multiple R-squared: 0.8205, Adjusted R-squared: 0.7871
F-statistic: 24.61 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,] 0.0037382126 0.0074764252 0.9962618
[2,] 0.0013027696 0.0026055392 0.9986972
[3,] 0.0005770839 0.0011541677 0.9994229
[4,] 0.0001684256 0.0003368511 0.9998316
[5,] 0.0013934864 0.0027869729 0.9986065
[6,] 0.0019627390 0.0039254779 0.9980373
[7,] 0.0011026917 0.0022053834 0.9988973
[8,] 0.0009424959 0.0018849918 0.9990575
[9,] 0.0032067035 0.0064134070 0.9967933
[10,] 0.0024060627 0.0048121253 0.9975939
[11,] 0.0012400359 0.0024800718 0.9987600
[12,] 0.0018628549 0.0037257099 0.9981371
[13,] 0.0010687328 0.0021374657 0.9989313
[14,] 0.0010157072 0.0020314144 0.9989843
[15,] 0.0015826985 0.0031653971 0.9984173
[16,] 0.0011117385 0.0022234771 0.9988883
[17,] 0.0022824223 0.0045648447 0.9977176
[18,] 0.0025828167 0.0051656335 0.9974172
[19,] 0.0021762141 0.0043524282 0.9978238
[20,] 0.0041591608 0.0083183215 0.9958408
[21,] 0.0039187133 0.0078374265 0.9960813
[22,] 0.0041918719 0.0083837439 0.9958081
[23,] 0.0034448014 0.0068896028 0.9965552
[24,] 0.0060131578 0.0120263157 0.9939868
[25,] 0.0081690517 0.0163381035 0.9918309
[26,] 0.0093671133 0.0187342266 0.9906329
[27,] 0.0219144196 0.0438288391 0.9780856
[28,] 0.0247620852 0.0495241705 0.9752379
[29,] 0.0291592616 0.0583185233 0.9708407
[30,] 0.0309648295 0.0619296590 0.9690352
[31,] 0.0285134837 0.0570269673 0.9714865
[32,] 0.0308305462 0.0616610924 0.9691695
[33,] 0.0299177394 0.0598354787 0.9700823
[34,] 0.0353197690 0.0706395380 0.9646802
[35,] 0.0422650960 0.0845301920 0.9577349
[36,] 0.0650175434 0.1300350868 0.9349825
[37,] 0.1365364062 0.2730728124 0.8634636
[38,] 0.3061731520 0.6123463040 0.6938268
[39,] 0.5425172768 0.9149654464 0.4574827
[40,] 0.6271950817 0.7456098366 0.3728049
[41,] 0.8507707270 0.2984585459 0.1492293
[42,] 0.8612456085 0.2775087830 0.1387544
[43,] 0.8749474989 0.2501050022 0.1250525
[44,] 0.8316966835 0.3366066331 0.1683033
[45,] 0.7968545140 0.4062909720 0.2031455
[46,] 0.7087164441 0.5825671117 0.2912836
[47,] 0.6310561195 0.7378877610 0.3689439
[48,] 0.6064909266 0.7870181468 0.3935091
[49,] 0.8805204346 0.2389591308 0.1194796
[50,] 0.8159654723 0.3680690554 0.1840345
[51,] 0.8165711420 0.3668577160 0.1834289
> postscript(file="/var/www/html/rcomp/tmp/15vs91229524386.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/20bpc1229524386.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/3dzq31229524386.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/4f7d91229524386.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/58h451229524386.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 7
14176.374 12330.111 10405.024 7443.058 9968.504 2185.640 -1866.986
8 9 10 11 12 13 14
5998.178 -5606.375 -4567.521 -2919.451 -6111.410 -11592.057 -12442.928
15 16 17 18 19 20 21
-10596.606 -11278.152 -11031.625 -13970.134 -17115.618 -15445.889 -7810.879
22 23 24 25 26 27 28
-3100.178 -2518.886 -1980.350 -6805.420 -8014.577 -6934.900 -7646.005
29 30 31 32 33 34 35
-8227.357 -9997.026 -4750.230 -7550.597 -5051.444 -2652.819 2235.517
36 37 38 39 40 41 42
3880.537 -2096.070 -2818.339 -3911.248 1759.295 1204.348 1771.637
43 44 45 46 47 48 49
11053.763 8320.232 7746.891 7186.548 4286.364 11071.603 4345.078
50 51 52 53 54 55 56
4721.009 7209.046 5193.439 3886.598 9856.450 11424.706 17757.791
57 58 59 60 61 62 63
15811.087 23935.678 22498.385 25428.602 12909.375 16453.215 14562.265
64 65 66 67 68 69 70
12411.103 9726.532 18182.099 14707.441 15281.151 13056.068 7810.353
71 72 73 74 75 76 77
8646.886 7255.631 -10937.281 -10228.490 -10733.582 -7882.738 -5527.000
78 79 80 81 82 83 84
-8028.665 -13453.075 -24360.865 -18145.348 -28612.061 -32228.815 -39544.613
> postscript(file="/var/www/html/rcomp/tmp/6pirx1229524386.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 14176.374 NA
1 12330.111 14176.374
2 10405.024 12330.111
3 7443.058 10405.024
4 9968.504 7443.058
5 2185.640 9968.504
6 -1866.986 2185.640
7 5998.178 -1866.986
8 -5606.375 5998.178
9 -4567.521 -5606.375
10 -2919.451 -4567.521
11 -6111.410 -2919.451
12 -11592.057 -6111.410
13 -12442.928 -11592.057
14 -10596.606 -12442.928
15 -11278.152 -10596.606
16 -11031.625 -11278.152
17 -13970.134 -11031.625
18 -17115.618 -13970.134
19 -15445.889 -17115.618
20 -7810.879 -15445.889
21 -3100.178 -7810.879
22 -2518.886 -3100.178
23 -1980.350 -2518.886
24 -6805.420 -1980.350
25 -8014.577 -6805.420
26 -6934.900 -8014.577
27 -7646.005 -6934.900
28 -8227.357 -7646.005
29 -9997.026 -8227.357
30 -4750.230 -9997.026
31 -7550.597 -4750.230
32 -5051.444 -7550.597
33 -2652.819 -5051.444
34 2235.517 -2652.819
35 3880.537 2235.517
36 -2096.070 3880.537
37 -2818.339 -2096.070
38 -3911.248 -2818.339
39 1759.295 -3911.248
40 1204.348 1759.295
41 1771.637 1204.348
42 11053.763 1771.637
43 8320.232 11053.763
44 7746.891 8320.232
45 7186.548 7746.891
46 4286.364 7186.548
47 11071.603 4286.364
48 4345.078 11071.603
49 4721.009 4345.078
50 7209.046 4721.009
51 5193.439 7209.046
52 3886.598 5193.439
53 9856.450 3886.598
54 11424.706 9856.450
55 17757.791 11424.706
56 15811.087 17757.791
57 23935.678 15811.087
58 22498.385 23935.678
59 25428.602 22498.385
60 12909.375 25428.602
61 16453.215 12909.375
62 14562.265 16453.215
63 12411.103 14562.265
64 9726.532 12411.103
65 18182.099 9726.532
66 14707.441 18182.099
67 15281.151 14707.441
68 13056.068 15281.151
69 7810.353 13056.068
70 8646.886 7810.353
71 7255.631 8646.886
72 -10937.281 7255.631
73 -10228.490 -10937.281
74 -10733.582 -10228.490
75 -7882.738 -10733.582
76 -5527.000 -7882.738
77 -8028.665 -5527.000
78 -13453.075 -8028.665
79 -24360.865 -13453.075
80 -18145.348 -24360.865
81 -28612.061 -18145.348
82 -32228.815 -28612.061
83 -39544.613 -32228.815
84 NA -39544.613
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 12330.111 14176.374
[2,] 10405.024 12330.111
[3,] 7443.058 10405.024
[4,] 9968.504 7443.058
[5,] 2185.640 9968.504
[6,] -1866.986 2185.640
[7,] 5998.178 -1866.986
[8,] -5606.375 5998.178
[9,] -4567.521 -5606.375
[10,] -2919.451 -4567.521
[11,] -6111.410 -2919.451
[12,] -11592.057 -6111.410
[13,] -12442.928 -11592.057
[14,] -10596.606 -12442.928
[15,] -11278.152 -10596.606
[16,] -11031.625 -11278.152
[17,] -13970.134 -11031.625
[18,] -17115.618 -13970.134
[19,] -15445.889 -17115.618
[20,] -7810.879 -15445.889
[21,] -3100.178 -7810.879
[22,] -2518.886 -3100.178
[23,] -1980.350 -2518.886
[24,] -6805.420 -1980.350
[25,] -8014.577 -6805.420
[26,] -6934.900 -8014.577
[27,] -7646.005 -6934.900
[28,] -8227.357 -7646.005
[29,] -9997.026 -8227.357
[30,] -4750.230 -9997.026
[31,] -7550.597 -4750.230
[32,] -5051.444 -7550.597
[33,] -2652.819 -5051.444
[34,] 2235.517 -2652.819
[35,] 3880.537 2235.517
[36,] -2096.070 3880.537
[37,] -2818.339 -2096.070
[38,] -3911.248 -2818.339
[39,] 1759.295 -3911.248
[40,] 1204.348 1759.295
[41,] 1771.637 1204.348
[42,] 11053.763 1771.637
[43,] 8320.232 11053.763
[44,] 7746.891 8320.232
[45,] 7186.548 7746.891
[46,] 4286.364 7186.548
[47,] 11071.603 4286.364
[48,] 4345.078 11071.603
[49,] 4721.009 4345.078
[50,] 7209.046 4721.009
[51,] 5193.439 7209.046
[52,] 3886.598 5193.439
[53,] 9856.450 3886.598
[54,] 11424.706 9856.450
[55,] 17757.791 11424.706
[56,] 15811.087 17757.791
[57,] 23935.678 15811.087
[58,] 22498.385 23935.678
[59,] 25428.602 22498.385
[60,] 12909.375 25428.602
[61,] 16453.215 12909.375
[62,] 14562.265 16453.215
[63,] 12411.103 14562.265
[64,] 9726.532 12411.103
[65,] 18182.099 9726.532
[66,] 14707.441 18182.099
[67,] 15281.151 14707.441
[68,] 13056.068 15281.151
[69,] 7810.353 13056.068
[70,] 8646.886 7810.353
[71,] 7255.631 8646.886
[72,] -10937.281 7255.631
[73,] -10228.490 -10937.281
[74,] -10733.582 -10228.490
[75,] -7882.738 -10733.582
[76,] -5527.000 -7882.738
[77,] -8028.665 -5527.000
[78,] -13453.075 -8028.665
[79,] -24360.865 -13453.075
[80,] -18145.348 -24360.865
[81,] -28612.061 -18145.348
[82,] -32228.815 -28612.061
[83,] -39544.613 -32228.815
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 12330.111 14176.374
2 10405.024 12330.111
3 7443.058 10405.024
4 9968.504 7443.058
5 2185.640 9968.504
6 -1866.986 2185.640
7 5998.178 -1866.986
8 -5606.375 5998.178
9 -4567.521 -5606.375
10 -2919.451 -4567.521
11 -6111.410 -2919.451
12 -11592.057 -6111.410
13 -12442.928 -11592.057
14 -10596.606 -12442.928
15 -11278.152 -10596.606
16 -11031.625 -11278.152
17 -13970.134 -11031.625
18 -17115.618 -13970.134
19 -15445.889 -17115.618
20 -7810.879 -15445.889
21 -3100.178 -7810.879
22 -2518.886 -3100.178
23 -1980.350 -2518.886
24 -6805.420 -1980.350
25 -8014.577 -6805.420
26 -6934.900 -8014.577
27 -7646.005 -6934.900
28 -8227.357 -7646.005
29 -9997.026 -8227.357
30 -4750.230 -9997.026
31 -7550.597 -4750.230
32 -5051.444 -7550.597
33 -2652.819 -5051.444
34 2235.517 -2652.819
35 3880.537 2235.517
36 -2096.070 3880.537
37 -2818.339 -2096.070
38 -3911.248 -2818.339
39 1759.295 -3911.248
40 1204.348 1759.295
41 1771.637 1204.348
42 11053.763 1771.637
43 8320.232 11053.763
44 7746.891 8320.232
45 7186.548 7746.891
46 4286.364 7186.548
47 11071.603 4286.364
48 4345.078 11071.603
49 4721.009 4345.078
50 7209.046 4721.009
51 5193.439 7209.046
52 3886.598 5193.439
53 9856.450 3886.598
54 11424.706 9856.450
55 17757.791 11424.706
56 15811.087 17757.791
57 23935.678 15811.087
58 22498.385 23935.678
59 25428.602 22498.385
60 12909.375 25428.602
61 16453.215 12909.375
62 14562.265 16453.215
63 12411.103 14562.265
64 9726.532 12411.103
65 18182.099 9726.532
66 14707.441 18182.099
67 15281.151 14707.441
68 13056.068 15281.151
69 7810.353 13056.068
70 8646.886 7810.353
71 7255.631 8646.886
72 -10937.281 7255.631
73 -10228.490 -10937.281
74 -10733.582 -10228.490
75 -7882.738 -10733.582
76 -5527.000 -7882.738
77 -8028.665 -5527.000
78 -13453.075 -8028.665
79 -24360.865 -13453.075
80 -18145.348 -24360.865
81 -28612.061 -18145.348
82 -32228.815 -28612.061
83 -39544.613 -32228.815
> 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/7hcwk1229524386.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/8h1io1229524386.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/9fxm11229524386.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/10dmpe1229524386.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/11whx31229524386.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/120msi1229524386.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/1399wf1229524386.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/1421h51229524386.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/15jio41229524386.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/16l6ui1229524386.tab")
+ }
>
> system("convert tmp/15vs91229524386.ps tmp/15vs91229524386.png")
> system("convert tmp/20bpc1229524386.ps tmp/20bpc1229524386.png")
> system("convert tmp/3dzq31229524386.ps tmp/3dzq31229524386.png")
> system("convert tmp/4f7d91229524386.ps tmp/4f7d91229524386.png")
> system("convert tmp/58h451229524386.ps tmp/58h451229524386.png")
> system("convert tmp/6pirx1229524386.ps tmp/6pirx1229524386.png")
> system("convert tmp/7hcwk1229524386.ps tmp/7hcwk1229524386.png")
> system("convert tmp/8h1io1229524386.ps tmp/8h1io1229524386.png")
> system("convert tmp/9fxm11229524386.ps tmp/9fxm11229524386.png")
> system("convert tmp/10dmpe1229524386.ps tmp/10dmpe1229524386.png")
>
>
> proc.time()
user system elapsed
5.576 2.790 5.953