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(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 = 'No Linear Trend'
> par2 = 'Do not include Seasonal 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
1 180144 1235.8
2 173666 1147.1
3 165688 1376.9
4 161570 1157.7
5 156145 1506.0
6 153730 1271.3
7 182698 1240.2
8 200765 1408.3
9 176512 1334.6
10 166618 1601.2
11 158644 1566.4
12 159585 1297.5
13 163095 1487.6
14 159044 1320.9
15 155511 1514.0
16 153745 1290.9
17 150569 1392.5
18 150605 1288.2
19 179612 1304.4
20 194690 1297.8
21 189917 1211.0
22 184128 1454.0
23 175335 1405.7
24 179566 1160.8
25 181140 1492.1
26 177876 1263.0
27 175041 1376.3
28 169292 1368.6
29 166070 1427.6
30 166972 1339.8
31 206348 1248.3
32 215706 1309.8
33 202108 1424.0
34 195411 1590.5
35 193111 1423.1
36 195198 1355.3
37 198770 1515.0
38 194163 1385.6
39 190420 1430.0
40 189733 1494.2
41 186029 1580.9
42 191531 1369.8
43 232571 1407.5
44 243477 1388.3
45 227247 1478.5
46 217859 1630.4
47 208679 1413.5
48 213188 1493.8
49 216234 1641.3
50 213586 1465.0
51 209465 1725.1
52 204045 1628.4
53 200237 1679.8
54 203666 1876.0
55 241476 1669.4
56 260307 1712.4
57 243324 1768.8
58 244460 1820.5
59 233575 1776.2
60 237217 1693.7
61 235243 1799.1
62 230354 1917.5
63 227184 1887.2
64 221678 1787.8
65 217142 1803.8
66 219452 2196.4
67 256446 1759.5
68 265845 2002.6
69 248624 2056.8
70 241114 1851.1
71 229245 1984.3
72 231805 1725.3
73 219277 2096.6
74 219313 1792.2
75 212610 2029.9
76 214771 1785.3
77 211142 2026.5
78 211457 1930.8
79 240048 1845.5
80 240636 1943.1
81 230580 2066.8
82 208795 2354.4
83 197922 2190.7
84 194596 1929.6
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Azi\353`
98514.99 64.83
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42610 -16723 1152 16310 54963
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 98514.995 14554.952 6.768 1.80e-09 ***
`Azi\353` 64.827 9.005 7.199 2.64e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23190 on 82 degrees of freedom
Multiple R-squared: 0.3872, Adjusted R-squared: 0.3798
F-statistic: 51.82 on 1 and 82 DF, p-value: 2.635e-10
> 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.07051265 0.14102530 0.92948735
[2,] 0.05762735 0.11525471 0.94237265
[3,] 0.05217776 0.10435552 0.94782224
[4,] 0.21378268 0.42756537 0.78621732
[5,] 0.13412533 0.26825066 0.86587467
[6,] 0.08751455 0.17502911 0.91248545
[7,] 0.06659068 0.13318136 0.93340932
[8,] 0.04969393 0.09938785 0.95030607
[9,] 0.03236713 0.06473426 0.96763287
[10,] 0.02419105 0.04838209 0.97580895
[11,] 0.02086988 0.04173976 0.97913012
[12,] 0.02148670 0.04297341 0.97851330
[13,] 0.02619870 0.05239740 0.97380130
[14,] 0.03286705 0.06573409 0.96713295
[15,] 0.02813821 0.05627642 0.97186179
[16,] 0.04987125 0.09974249 0.95012875
[17,] 0.04940654 0.09881309 0.95059346
[18,] 0.05143674 0.10287348 0.94856326
[19,] 0.04173096 0.08346191 0.95826904
[20,] 0.02922830 0.05845659 0.97077170
[21,] 0.02848960 0.05697920 0.97151040
[22,] 0.02122351 0.04244702 0.97877649
[23,] 0.01746028 0.03492056 0.98253972
[24,] 0.01603546 0.03207093 0.98396454
[25,] 0.01827799 0.03655598 0.98172201
[26,] 0.02157808 0.04315615 0.97842192
[27,] 0.05005785 0.10011569 0.94994215
[28,] 0.14652479 0.29304958 0.85347521
[29,] 0.19621279 0.39242559 0.80378721
[30,] 0.23615710 0.47231421 0.76384290
[31,] 0.24315907 0.48631814 0.75684093
[32,] 0.24770656 0.49541311 0.75229344
[33,] 0.27067095 0.54134189 0.72932905
[34,] 0.27724266 0.55448532 0.72275734
[35,] 0.29365827 0.58731653 0.70634173
[36,] 0.32555895 0.65111790 0.67444105
[37,] 0.39109180 0.78218359 0.60890820
[38,] 0.46507760 0.93015521 0.53492240
[39,] 0.67495873 0.65008254 0.32504127
[40,] 0.87040141 0.25919717 0.12959859
[41,] 0.90040650 0.19918700 0.09959350
[42,] 0.89752493 0.20495014 0.10247507
[43,] 0.89985007 0.20029986 0.10014993
[44,] 0.89938400 0.20123200 0.10061600
[45,] 0.88947318 0.22105364 0.11052682
[46,] 0.89511956 0.20976087 0.10488044
[47,] 0.88991912 0.22016176 0.11008088
[48,] 0.91376036 0.17247928 0.08623964
[49,] 0.95101681 0.09796638 0.04898319
[50,] 0.96049458 0.07901083 0.03950542
[51,] 0.96162535 0.07674930 0.03837465
[52,] 0.98156100 0.03687801 0.01843900
[53,] 0.97838936 0.04322127 0.02161064
[54,] 0.97502891 0.04994219 0.02497109
[55,] 0.96355910 0.07288180 0.03644090
[56,] 0.95110513 0.09778973 0.04889487
[57,] 0.93109136 0.13781727 0.06890864
[58,] 0.90276992 0.19446016 0.09723008
[59,] 0.86507800 0.26984400 0.13492200
[60,] 0.82694272 0.34611456 0.17305728
[61,] 0.79761997 0.40476006 0.20238003
[62,] 0.76384866 0.47230268 0.23615134
[63,] 0.79936158 0.40127683 0.20063842
[64,] 0.94343390 0.11313219 0.05656610
[65,] 0.97435223 0.05129554 0.02564777
[66,] 0.97326192 0.05347617 0.02673808
[67,] 0.96198564 0.07602872 0.03801436
[68,] 0.93900718 0.12198565 0.06099282
[69,] 0.90811163 0.18377675 0.09188837
[70,] 0.85255114 0.29489773 0.14744886
[71,] 0.78251381 0.43497237 0.21748619
[72,] 0.70918247 0.58163506 0.29081753
[73,] 0.60393308 0.79213384 0.39606692
[74,] 0.50486594 0.99026813 0.49513406
[75,] 0.39449314 0.78898629 0.60550686
> postscript(file="/var/www/html/rcomp/tmp/10io81229522854.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/24ob71229522854.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/3btwc1229522854.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/4s06t1229522854.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/580x21229522854.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
1515.64733 787.81308 -22087.45963 -11995.35442 -39999.64112 -27199.71551
7 8 9 10 11 12
3784.40799 10953.96873 -8521.27235 -35698.18316 -41416.19930 -23043.18612
13 14 15 16 17 18
-31856.82207 -25101.14078 -41152.25809 -28455.32711 -38217.76274 -31420.29388
19 20 21 22 23 24
-3463.49326 12042.36575 12896.35996 -8645.63076 -14307.48075 5799.68150
25 26 27 28 29 30
-14103.54412 -2515.65040 -12695.56335 -17945.39451 -24992.19473 -18398.37339
31 32 33 34 35 36
26909.30830 32280.44028 11279.18291 -6211.53295 2340.52732 8822.80621
37 38 39 40 41 42
2041.91478 5823.54441 -797.77982 -5646.68107 -14971.19258 4215.81294
43 44 45 46 47 48
42811.83043 54962.51118 32885.10475 13649.86487 18530.86770 17834.24978
49 50 51 52 53 54
11318.24924 20099.27090 -883.26361 -34.48088 -7174.59497 -16464.67637
55 56 57 58 59 60
34738.60710 50782.04084 30142.79114 27927.22892 19914.07044 28904.30803
61 62 63 64 65 66
20097.52934 7532.99806 6327.25986 7265.07582 1691.84186 -21449.28635
67 68 69 70 71 72
43867.68338 37507.20995 16772.57992 22597.51898 2093.54629 21443.77096
73 74 75 76 77 78
-15154.53954 4614.83648 -17497.57049 520.14363 -18745.15827 -12226.20267
79 80 81 82 83 84
21894.55086 16155.42373 -1919.69130 -42348.97167 -42609.77175 -29009.41012
> postscript(file="/var/www/html/rcomp/tmp/6kyry1229522854.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 1515.64733 NA
1 787.81308 1515.64733
2 -22087.45963 787.81308
3 -11995.35442 -22087.45963
4 -39999.64112 -11995.35442
5 -27199.71551 -39999.64112
6 3784.40799 -27199.71551
7 10953.96873 3784.40799
8 -8521.27235 10953.96873
9 -35698.18316 -8521.27235
10 -41416.19930 -35698.18316
11 -23043.18612 -41416.19930
12 -31856.82207 -23043.18612
13 -25101.14078 -31856.82207
14 -41152.25809 -25101.14078
15 -28455.32711 -41152.25809
16 -38217.76274 -28455.32711
17 -31420.29388 -38217.76274
18 -3463.49326 -31420.29388
19 12042.36575 -3463.49326
20 12896.35996 12042.36575
21 -8645.63076 12896.35996
22 -14307.48075 -8645.63076
23 5799.68150 -14307.48075
24 -14103.54412 5799.68150
25 -2515.65040 -14103.54412
26 -12695.56335 -2515.65040
27 -17945.39451 -12695.56335
28 -24992.19473 -17945.39451
29 -18398.37339 -24992.19473
30 26909.30830 -18398.37339
31 32280.44028 26909.30830
32 11279.18291 32280.44028
33 -6211.53295 11279.18291
34 2340.52732 -6211.53295
35 8822.80621 2340.52732
36 2041.91478 8822.80621
37 5823.54441 2041.91478
38 -797.77982 5823.54441
39 -5646.68107 -797.77982
40 -14971.19258 -5646.68107
41 4215.81294 -14971.19258
42 42811.83043 4215.81294
43 54962.51118 42811.83043
44 32885.10475 54962.51118
45 13649.86487 32885.10475
46 18530.86770 13649.86487
47 17834.24978 18530.86770
48 11318.24924 17834.24978
49 20099.27090 11318.24924
50 -883.26361 20099.27090
51 -34.48088 -883.26361
52 -7174.59497 -34.48088
53 -16464.67637 -7174.59497
54 34738.60710 -16464.67637
55 50782.04084 34738.60710
56 30142.79114 50782.04084
57 27927.22892 30142.79114
58 19914.07044 27927.22892
59 28904.30803 19914.07044
60 20097.52934 28904.30803
61 7532.99806 20097.52934
62 6327.25986 7532.99806
63 7265.07582 6327.25986
64 1691.84186 7265.07582
65 -21449.28635 1691.84186
66 43867.68338 -21449.28635
67 37507.20995 43867.68338
68 16772.57992 37507.20995
69 22597.51898 16772.57992
70 2093.54629 22597.51898
71 21443.77096 2093.54629
72 -15154.53954 21443.77096
73 4614.83648 -15154.53954
74 -17497.57049 4614.83648
75 520.14363 -17497.57049
76 -18745.15827 520.14363
77 -12226.20267 -18745.15827
78 21894.55086 -12226.20267
79 16155.42373 21894.55086
80 -1919.69130 16155.42373
81 -42348.97167 -1919.69130
82 -42609.77175 -42348.97167
83 -29009.41012 -42609.77175
84 NA -29009.41012
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 787.81308 1515.64733
[2,] -22087.45963 787.81308
[3,] -11995.35442 -22087.45963
[4,] -39999.64112 -11995.35442
[5,] -27199.71551 -39999.64112
[6,] 3784.40799 -27199.71551
[7,] 10953.96873 3784.40799
[8,] -8521.27235 10953.96873
[9,] -35698.18316 -8521.27235
[10,] -41416.19930 -35698.18316
[11,] -23043.18612 -41416.19930
[12,] -31856.82207 -23043.18612
[13,] -25101.14078 -31856.82207
[14,] -41152.25809 -25101.14078
[15,] -28455.32711 -41152.25809
[16,] -38217.76274 -28455.32711
[17,] -31420.29388 -38217.76274
[18,] -3463.49326 -31420.29388
[19,] 12042.36575 -3463.49326
[20,] 12896.35996 12042.36575
[21,] -8645.63076 12896.35996
[22,] -14307.48075 -8645.63076
[23,] 5799.68150 -14307.48075
[24,] -14103.54412 5799.68150
[25,] -2515.65040 -14103.54412
[26,] -12695.56335 -2515.65040
[27,] -17945.39451 -12695.56335
[28,] -24992.19473 -17945.39451
[29,] -18398.37339 -24992.19473
[30,] 26909.30830 -18398.37339
[31,] 32280.44028 26909.30830
[32,] 11279.18291 32280.44028
[33,] -6211.53295 11279.18291
[34,] 2340.52732 -6211.53295
[35,] 8822.80621 2340.52732
[36,] 2041.91478 8822.80621
[37,] 5823.54441 2041.91478
[38,] -797.77982 5823.54441
[39,] -5646.68107 -797.77982
[40,] -14971.19258 -5646.68107
[41,] 4215.81294 -14971.19258
[42,] 42811.83043 4215.81294
[43,] 54962.51118 42811.83043
[44,] 32885.10475 54962.51118
[45,] 13649.86487 32885.10475
[46,] 18530.86770 13649.86487
[47,] 17834.24978 18530.86770
[48,] 11318.24924 17834.24978
[49,] 20099.27090 11318.24924
[50,] -883.26361 20099.27090
[51,] -34.48088 -883.26361
[52,] -7174.59497 -34.48088
[53,] -16464.67637 -7174.59497
[54,] 34738.60710 -16464.67637
[55,] 50782.04084 34738.60710
[56,] 30142.79114 50782.04084
[57,] 27927.22892 30142.79114
[58,] 19914.07044 27927.22892
[59,] 28904.30803 19914.07044
[60,] 20097.52934 28904.30803
[61,] 7532.99806 20097.52934
[62,] 6327.25986 7532.99806
[63,] 7265.07582 6327.25986
[64,] 1691.84186 7265.07582
[65,] -21449.28635 1691.84186
[66,] 43867.68338 -21449.28635
[67,] 37507.20995 43867.68338
[68,] 16772.57992 37507.20995
[69,] 22597.51898 16772.57992
[70,] 2093.54629 22597.51898
[71,] 21443.77096 2093.54629
[72,] -15154.53954 21443.77096
[73,] 4614.83648 -15154.53954
[74,] -17497.57049 4614.83648
[75,] 520.14363 -17497.57049
[76,] -18745.15827 520.14363
[77,] -12226.20267 -18745.15827
[78,] 21894.55086 -12226.20267
[79,] 16155.42373 21894.55086
[80,] -1919.69130 16155.42373
[81,] -42348.97167 -1919.69130
[82,] -42609.77175 -42348.97167
[83,] -29009.41012 -42609.77175
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 787.81308 1515.64733
2 -22087.45963 787.81308
3 -11995.35442 -22087.45963
4 -39999.64112 -11995.35442
5 -27199.71551 -39999.64112
6 3784.40799 -27199.71551
7 10953.96873 3784.40799
8 -8521.27235 10953.96873
9 -35698.18316 -8521.27235
10 -41416.19930 -35698.18316
11 -23043.18612 -41416.19930
12 -31856.82207 -23043.18612
13 -25101.14078 -31856.82207
14 -41152.25809 -25101.14078
15 -28455.32711 -41152.25809
16 -38217.76274 -28455.32711
17 -31420.29388 -38217.76274
18 -3463.49326 -31420.29388
19 12042.36575 -3463.49326
20 12896.35996 12042.36575
21 -8645.63076 12896.35996
22 -14307.48075 -8645.63076
23 5799.68150 -14307.48075
24 -14103.54412 5799.68150
25 -2515.65040 -14103.54412
26 -12695.56335 -2515.65040
27 -17945.39451 -12695.56335
28 -24992.19473 -17945.39451
29 -18398.37339 -24992.19473
30 26909.30830 -18398.37339
31 32280.44028 26909.30830
32 11279.18291 32280.44028
33 -6211.53295 11279.18291
34 2340.52732 -6211.53295
35 8822.80621 2340.52732
36 2041.91478 8822.80621
37 5823.54441 2041.91478
38 -797.77982 5823.54441
39 -5646.68107 -797.77982
40 -14971.19258 -5646.68107
41 4215.81294 -14971.19258
42 42811.83043 4215.81294
43 54962.51118 42811.83043
44 32885.10475 54962.51118
45 13649.86487 32885.10475
46 18530.86770 13649.86487
47 17834.24978 18530.86770
48 11318.24924 17834.24978
49 20099.27090 11318.24924
50 -883.26361 20099.27090
51 -34.48088 -883.26361
52 -7174.59497 -34.48088
53 -16464.67637 -7174.59497
54 34738.60710 -16464.67637
55 50782.04084 34738.60710
56 30142.79114 50782.04084
57 27927.22892 30142.79114
58 19914.07044 27927.22892
59 28904.30803 19914.07044
60 20097.52934 28904.30803
61 7532.99806 20097.52934
62 6327.25986 7532.99806
63 7265.07582 6327.25986
64 1691.84186 7265.07582
65 -21449.28635 1691.84186
66 43867.68338 -21449.28635
67 37507.20995 43867.68338
68 16772.57992 37507.20995
69 22597.51898 16772.57992
70 2093.54629 22597.51898
71 21443.77096 2093.54629
72 -15154.53954 21443.77096
73 4614.83648 -15154.53954
74 -17497.57049 4614.83648
75 520.14363 -17497.57049
76 -18745.15827 520.14363
77 -12226.20267 -18745.15827
78 21894.55086 -12226.20267
79 16155.42373 21894.55086
80 -1919.69130 16155.42373
81 -42348.97167 -1919.69130
82 -42609.77175 -42348.97167
83 -29009.41012 -42609.77175
> 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/7zonw1229522854.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/8d0mi1229522854.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/96a8s1229522854.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/10i0c41229522854.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/11dqli1229522854.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/12dn1j1229522854.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/13j25r1229522854.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/14oikc1229522855.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/15fxlh1229522855.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/16j00m1229522855.tab")
+ }
>
> system("convert tmp/10io81229522854.ps tmp/10io81229522854.png")
> system("convert tmp/24ob71229522854.ps tmp/24ob71229522854.png")
> system("convert tmp/3btwc1229522854.ps tmp/3btwc1229522854.png")
> system("convert tmp/4s06t1229522854.ps tmp/4s06t1229522854.png")
> system("convert tmp/580x21229522854.ps tmp/580x21229522854.png")
> system("convert tmp/6kyry1229522854.ps tmp/6kyry1229522854.png")
> system("convert tmp/7zonw1229522854.ps tmp/7zonw1229522854.png")
> system("convert tmp/8d0mi1229522854.ps tmp/8d0mi1229522854.png")
> system("convert tmp/96a8s1229522854.ps tmp/96a8s1229522854.png")
> system("convert tmp/10i0c41229522854.ps tmp/10i0c41229522854.png")
>
>
> proc.time()
user system elapsed
3.074 1.763 4.885