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 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
werkloosheid Azi\353 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 180144 1235.8 1 0 0 0 0 0 0 0 0 0 0
2 173666 1147.1 0 1 0 0 0 0 0 0 0 0 0
3 165688 1376.9 0 0 1 0 0 0 0 0 0 0 0
4 161570 1157.7 0 0 0 1 0 0 0 0 0 0 0
5 156145 1506.0 0 0 0 0 1 0 0 0 0 0 0
6 153730 1271.3 0 0 0 0 0 1 0 0 0 0 0
7 182698 1240.2 0 0 0 0 0 0 1 0 0 0 0
8 200765 1408.3 0 0 0 0 0 0 0 1 0 0 0
9 176512 1334.6 0 0 0 0 0 0 0 0 1 0 0
10 166618 1601.2 0 0 0 0 0 0 0 0 0 1 0
11 158644 1566.4 0 0 0 0 0 0 0 0 0 0 1
12 159585 1297.5 0 0 0 0 0 0 0 0 0 0 0
13 163095 1487.6 1 0 0 0 0 0 0 0 0 0 0
14 159044 1320.9 0 1 0 0 0 0 0 0 0 0 0
15 155511 1514.0 0 0 1 0 0 0 0 0 0 0 0
16 153745 1290.9 0 0 0 1 0 0 0 0 0 0 0
17 150569 1392.5 0 0 0 0 1 0 0 0 0 0 0
18 150605 1288.2 0 0 0 0 0 1 0 0 0 0 0
19 179612 1304.4 0 0 0 0 0 0 1 0 0 0 0
20 194690 1297.8 0 0 0 0 0 0 0 1 0 0 0
21 189917 1211.0 0 0 0 0 0 0 0 0 1 0 0
22 184128 1454.0 0 0 0 0 0 0 0 0 0 1 0
23 175335 1405.7 0 0 0 0 0 0 0 0 0 0 1
24 179566 1160.8 0 0 0 0 0 0 0 0 0 0 0
25 181140 1492.1 1 0 0 0 0 0 0 0 0 0 0
26 177876 1263.0 0 1 0 0 0 0 0 0 0 0 0
27 175041 1376.3 0 0 1 0 0 0 0 0 0 0 0
28 169292 1368.6 0 0 0 1 0 0 0 0 0 0 0
29 166070 1427.6 0 0 0 0 1 0 0 0 0 0 0
30 166972 1339.8 0 0 0 0 0 1 0 0 0 0 0
31 206348 1248.3 0 0 0 0 0 0 1 0 0 0 0
32 215706 1309.8 0 0 0 0 0 0 0 1 0 0 0
33 202108 1424.0 0 0 0 0 0 0 0 0 1 0 0
34 195411 1590.5 0 0 0 0 0 0 0 0 0 1 0
35 193111 1423.1 0 0 0 0 0 0 0 0 0 0 1
36 195198 1355.3 0 0 0 0 0 0 0 0 0 0 0
37 198770 1515.0 1 0 0 0 0 0 0 0 0 0 0
38 194163 1385.6 0 1 0 0 0 0 0 0 0 0 0
39 190420 1430.0 0 0 1 0 0 0 0 0 0 0 0
40 189733 1494.2 0 0 0 1 0 0 0 0 0 0 0
41 186029 1580.9 0 0 0 0 1 0 0 0 0 0 0
42 191531 1369.8 0 0 0 0 0 1 0 0 0 0 0
43 232571 1407.5 0 0 0 0 0 0 1 0 0 0 0
44 243477 1388.3 0 0 0 0 0 0 0 1 0 0 0
45 227247 1478.5 0 0 0 0 0 0 0 0 1 0 0
46 217859 1630.4 0 0 0 0 0 0 0 0 0 1 0
47 208679 1413.5 0 0 0 0 0 0 0 0 0 0 1
48 213188 1493.8 0 0 0 0 0 0 0 0 0 0 0
49 216234 1641.3 1 0 0 0 0 0 0 0 0 0 0
50 213586 1465.0 0 1 0 0 0 0 0 0 0 0 0
51 209465 1725.1 0 0 1 0 0 0 0 0 0 0 0
52 204045 1628.4 0 0 0 1 0 0 0 0 0 0 0
53 200237 1679.8 0 0 0 0 1 0 0 0 0 0 0
54 203666 1876.0 0 0 0 0 0 1 0 0 0 0 0
55 241476 1669.4 0 0 0 0 0 0 1 0 0 0 0
56 260307 1712.4 0 0 0 0 0 0 0 1 0 0 0
57 243324 1768.8 0 0 0 0 0 0 0 0 1 0 0
58 244460 1820.5 0 0 0 0 0 0 0 0 0 1 0
59 233575 1776.2 0 0 0 0 0 0 0 0 0 0 1
60 237217 1693.7 0 0 0 0 0 0 0 0 0 0 0
61 235243 1799.1 1 0 0 0 0 0 0 0 0 0 0
62 230354 1917.5 0 1 0 0 0 0 0 0 0 0 0
63 227184 1887.2 0 0 1 0 0 0 0 0 0 0 0
64 221678 1787.8 0 0 0 1 0 0 0 0 0 0 0
65 217142 1803.8 0 0 0 0 1 0 0 0 0 0 0
66 219452 2196.4 0 0 0 0 0 1 0 0 0 0 0
67 256446 1759.5 0 0 0 0 0 0 1 0 0 0 0
68 265845 2002.6 0 0 0 0 0 0 0 1 0 0 0
69 248624 2056.8 0 0 0 0 0 0 0 0 1 0 0
70 241114 1851.1 0 0 0 0 0 0 0 0 0 1 0
71 229245 1984.3 0 0 0 0 0 0 0 0 0 0 1
72 231805 1725.3 0 0 0 0 0 0 0 0 0 0 0
73 219277 2096.6 1 0 0 0 0 0 0 0 0 0 0
74 219313 1792.2 0 1 0 0 0 0 0 0 0 0 0
75 212610 2029.9 0 0 1 0 0 0 0 0 0 0 0
76 214771 1785.3 0 0 0 1 0 0 0 0 0 0 0
77 211142 2026.5 0 0 0 0 1 0 0 0 0 0 0
78 211457 1930.8 0 0 0 0 0 1 0 0 0 0 0
79 240048 1845.5 0 0 0 0 0 0 1 0 0 0 0
80 240636 1943.1 0 0 0 0 0 0 0 1 0 0 0
81 230580 2066.8 0 0 0 0 0 0 0 0 1 0 0
82 208795 2354.4 0 0 0 0 0 0 0 0 0 1 0
83 197922 2190.7 0 0 0 0 0 0 0 0 0 0 1
84 194596 1929.6 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Azi\353` M1 M2 M3 M4
94175.95 70.56 -8628.79 -2488.36 -17637.00 -12317.63
M5 M6 M7 M8 M9 M10
-25360.97 -22461.46 20118.58 25943.01 8408.05 -9846.34
M11
-13219.84
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-41669 -10273 1633 13133 31670
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 94175.954 13424.102 7.015 1.10e-09 ***
`Azi\353` 70.563 7.517 9.387 4.54e-14 ***
M1 -8628.786 9947.342 -0.867 0.3886
M2 -2488.362 9933.367 -0.251 0.8029
M3 -17637.001 9952.738 -1.772 0.0807 .
M4 -12317.626 9926.833 -1.241 0.2187
M5 -25360.969 9959.238 -2.546 0.0131 *
M6 -22461.458 9947.684 -2.258 0.0270 *
M7 20118.583 9927.551 2.027 0.0465 *
M8 25943.014 9935.229 2.611 0.0110 *
M9 8408.053 9952.825 0.845 0.4011
M10 -9846.339 10081.824 -0.977 0.3321
M11 -13219.844 9996.184 -1.322 0.1902
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 18570 on 71 degrees of freedom
Multiple R-squared: 0.6599, Adjusted R-squared: 0.6025
F-statistic: 11.48 on 12 and 71 DF, p-value: 2.215e-12
> 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.0005179284 0.0010358568 0.99948207
[2,] 0.0121404969 0.0242809938 0.98785950
[3,] 0.0032670941 0.0065341882 0.99673291
[4,] 0.0009755930 0.0019511860 0.99902441
[5,] 0.0018399950 0.0036799899 0.99816001
[6,] 0.0008538086 0.0017076172 0.99914619
[7,] 0.0005523896 0.0011047792 0.99944761
[8,] 0.0002646456 0.0005292912 0.99973535
[9,] 0.0002134976 0.0004269953 0.99978650
[10,] 0.0014576058 0.0029152115 0.99854239
[11,] 0.0019414725 0.0038829450 0.99805853
[12,] 0.0016618609 0.0033237217 0.99833814
[13,] 0.0055096956 0.0110193911 0.99449030
[14,] 0.0058876422 0.0117752844 0.99411236
[15,] 0.0103652606 0.0207305212 0.98963474
[16,] 0.0299353557 0.0598707115 0.97006464
[17,] 0.0384242145 0.0768484289 0.96157579
[18,] 0.0880115780 0.1760231561 0.91198842
[19,] 0.1417864931 0.2835729862 0.85821351
[20,] 0.1934804067 0.3869608134 0.80651959
[21,] 0.2980977836 0.5961955672 0.70190222
[22,] 0.3908756044 0.7817512088 0.60912440
[23,] 0.4825943789 0.9651887577 0.51740562
[24,] 0.5645269045 0.8709461910 0.43547310
[25,] 0.6488543088 0.7022913824 0.35114569
[26,] 0.7190668224 0.5618663552 0.28093318
[27,] 0.7941207411 0.4117585179 0.20587926
[28,] 0.8563487390 0.2873025221 0.14365126
[29,] 0.8963354831 0.2073290338 0.10366452
[30,] 0.9088215525 0.1823568950 0.09117845
[31,] 0.9241692411 0.1516615178 0.07583076
[32,] 0.9496274728 0.1007450545 0.05037253
[33,] 0.9434821701 0.1130356598 0.05651783
[34,] 0.9428838151 0.1142323698 0.05711618
[35,] 0.9532616855 0.0934766290 0.04673831
[36,] 0.9508681231 0.0982637538 0.04913188
[37,] 0.9519217263 0.0961565475 0.04807827
[38,] 0.9655283809 0.0689432382 0.03447162
[39,] 0.9731481373 0.0537037253 0.02685186
[40,] 0.9671394021 0.0657211958 0.03286060
[41,] 0.9523354686 0.0953290627 0.04766453
[42,] 0.9425807748 0.1148384505 0.05741923
[43,] 0.9234288120 0.1531423761 0.07657119
[44,] 0.8862027225 0.2275945549 0.11379728
[45,] 0.8728940112 0.2542119776 0.12710599
[46,] 0.8136136341 0.3727727318 0.18638637
[47,] 0.8036590509 0.3926818982 0.19634095
[48,] 0.7212261214 0.5575477571 0.27877388
[49,] 0.6237559415 0.7524881169 0.37624406
[50,] 0.5235421810 0.9529156381 0.47645782
[51,] 0.6216471674 0.7567056651 0.37835283
[52,] 0.4948086295 0.9896172591 0.50519137
[53,] 0.6231720955 0.7536558090 0.37682790
> postscript(file="/var/www/html/rcomp/tmp/1fc841229523647.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/2da4v1229523647.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/33e3h1229523647.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/418c31229523647.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/5bp3j1229523647.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
7394.6137 1035.1612 -8009.6633 -1979.5471 -18938.4272 -7691.7141
7 8 9 10 11 12
-19109.2348 -18728.3691 -20245.8869 -30697.6910 -32842.5801 -26146.9328
13 14 15 16 17 18
-27422.2441 -25850.7533 -27860.9020 -19203.5887 -16505.4842 -12009.2351
19 20 21 22 23 24
-26725.4035 -17006.1162 1880.7463 -2800.7622 -4812.0458 3480.0805
25 26 27 28 29 30
-9694.7793 -2933.1339 1385.6747 -9139.3629 -3481.2586 716.6947
31 32 33 34 35 36
3969.2019 3163.1233 -958.2526 -1149.6629 11736.1515 5387.5041
37 38 39 40 41 42
6319.3194 4702.7963 12975.4215 2438.8772 5660.3760 23158.7935
43 44 45 46 47 48
18958.5126 25394.8984 20335.0435 18482.8585 27981.5599 13604.4767
49 50 51 52 53 54
14871.1652 18523.0644 11197.1696 7281.2723 12889.6582 -425.3869
55 56 57 58 59 60
9382.9647 19355.3086 15927.4958 31669.7609 27284.2239 23527.8581
61 62 63 64 65 66
22745.2646 3361.1372 17477.8465 13666.4704 21044.7997 -7247.8922
67 68 69 70 71 72
17995.2046 4415.8172 905.2438 26164.5216 8269.9856 15886.0554
73 74 75 76 77 78
-14213.3394 1161.7281 -7165.5471 6935.8788 -669.6639 3498.7402
79 80 81 82 83 84
-4471.2456 -16594.6620 -17844.3899 -41669.0249 -37617.2950 -35739.0421
> postscript(file="/var/www/html/rcomp/tmp/65r8z1229523647.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 7394.6137 NA
1 1035.1612 7394.6137
2 -8009.6633 1035.1612
3 -1979.5471 -8009.6633
4 -18938.4272 -1979.5471
5 -7691.7141 -18938.4272
6 -19109.2348 -7691.7141
7 -18728.3691 -19109.2348
8 -20245.8869 -18728.3691
9 -30697.6910 -20245.8869
10 -32842.5801 -30697.6910
11 -26146.9328 -32842.5801
12 -27422.2441 -26146.9328
13 -25850.7533 -27422.2441
14 -27860.9020 -25850.7533
15 -19203.5887 -27860.9020
16 -16505.4842 -19203.5887
17 -12009.2351 -16505.4842
18 -26725.4035 -12009.2351
19 -17006.1162 -26725.4035
20 1880.7463 -17006.1162
21 -2800.7622 1880.7463
22 -4812.0458 -2800.7622
23 3480.0805 -4812.0458
24 -9694.7793 3480.0805
25 -2933.1339 -9694.7793
26 1385.6747 -2933.1339
27 -9139.3629 1385.6747
28 -3481.2586 -9139.3629
29 716.6947 -3481.2586
30 3969.2019 716.6947
31 3163.1233 3969.2019
32 -958.2526 3163.1233
33 -1149.6629 -958.2526
34 11736.1515 -1149.6629
35 5387.5041 11736.1515
36 6319.3194 5387.5041
37 4702.7963 6319.3194
38 12975.4215 4702.7963
39 2438.8772 12975.4215
40 5660.3760 2438.8772
41 23158.7935 5660.3760
42 18958.5126 23158.7935
43 25394.8984 18958.5126
44 20335.0435 25394.8984
45 18482.8585 20335.0435
46 27981.5599 18482.8585
47 13604.4767 27981.5599
48 14871.1652 13604.4767
49 18523.0644 14871.1652
50 11197.1696 18523.0644
51 7281.2723 11197.1696
52 12889.6582 7281.2723
53 -425.3869 12889.6582
54 9382.9647 -425.3869
55 19355.3086 9382.9647
56 15927.4958 19355.3086
57 31669.7609 15927.4958
58 27284.2239 31669.7609
59 23527.8581 27284.2239
60 22745.2646 23527.8581
61 3361.1372 22745.2646
62 17477.8465 3361.1372
63 13666.4704 17477.8465
64 21044.7997 13666.4704
65 -7247.8922 21044.7997
66 17995.2046 -7247.8922
67 4415.8172 17995.2046
68 905.2438 4415.8172
69 26164.5216 905.2438
70 8269.9856 26164.5216
71 15886.0554 8269.9856
72 -14213.3394 15886.0554
73 1161.7281 -14213.3394
74 -7165.5471 1161.7281
75 6935.8788 -7165.5471
76 -669.6639 6935.8788
77 3498.7402 -669.6639
78 -4471.2456 3498.7402
79 -16594.6620 -4471.2456
80 -17844.3899 -16594.6620
81 -41669.0249 -17844.3899
82 -37617.2950 -41669.0249
83 -35739.0421 -37617.2950
84 NA -35739.0421
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1035.1612 7394.6137
[2,] -8009.6633 1035.1612
[3,] -1979.5471 -8009.6633
[4,] -18938.4272 -1979.5471
[5,] -7691.7141 -18938.4272
[6,] -19109.2348 -7691.7141
[7,] -18728.3691 -19109.2348
[8,] -20245.8869 -18728.3691
[9,] -30697.6910 -20245.8869
[10,] -32842.5801 -30697.6910
[11,] -26146.9328 -32842.5801
[12,] -27422.2441 -26146.9328
[13,] -25850.7533 -27422.2441
[14,] -27860.9020 -25850.7533
[15,] -19203.5887 -27860.9020
[16,] -16505.4842 -19203.5887
[17,] -12009.2351 -16505.4842
[18,] -26725.4035 -12009.2351
[19,] -17006.1162 -26725.4035
[20,] 1880.7463 -17006.1162
[21,] -2800.7622 1880.7463
[22,] -4812.0458 -2800.7622
[23,] 3480.0805 -4812.0458
[24,] -9694.7793 3480.0805
[25,] -2933.1339 -9694.7793
[26,] 1385.6747 -2933.1339
[27,] -9139.3629 1385.6747
[28,] -3481.2586 -9139.3629
[29,] 716.6947 -3481.2586
[30,] 3969.2019 716.6947
[31,] 3163.1233 3969.2019
[32,] -958.2526 3163.1233
[33,] -1149.6629 -958.2526
[34,] 11736.1515 -1149.6629
[35,] 5387.5041 11736.1515
[36,] 6319.3194 5387.5041
[37,] 4702.7963 6319.3194
[38,] 12975.4215 4702.7963
[39,] 2438.8772 12975.4215
[40,] 5660.3760 2438.8772
[41,] 23158.7935 5660.3760
[42,] 18958.5126 23158.7935
[43,] 25394.8984 18958.5126
[44,] 20335.0435 25394.8984
[45,] 18482.8585 20335.0435
[46,] 27981.5599 18482.8585
[47,] 13604.4767 27981.5599
[48,] 14871.1652 13604.4767
[49,] 18523.0644 14871.1652
[50,] 11197.1696 18523.0644
[51,] 7281.2723 11197.1696
[52,] 12889.6582 7281.2723
[53,] -425.3869 12889.6582
[54,] 9382.9647 -425.3869
[55,] 19355.3086 9382.9647
[56,] 15927.4958 19355.3086
[57,] 31669.7609 15927.4958
[58,] 27284.2239 31669.7609
[59,] 23527.8581 27284.2239
[60,] 22745.2646 23527.8581
[61,] 3361.1372 22745.2646
[62,] 17477.8465 3361.1372
[63,] 13666.4704 17477.8465
[64,] 21044.7997 13666.4704
[65,] -7247.8922 21044.7997
[66,] 17995.2046 -7247.8922
[67,] 4415.8172 17995.2046
[68,] 905.2438 4415.8172
[69,] 26164.5216 905.2438
[70,] 8269.9856 26164.5216
[71,] 15886.0554 8269.9856
[72,] -14213.3394 15886.0554
[73,] 1161.7281 -14213.3394
[74,] -7165.5471 1161.7281
[75,] 6935.8788 -7165.5471
[76,] -669.6639 6935.8788
[77,] 3498.7402 -669.6639
[78,] -4471.2456 3498.7402
[79,] -16594.6620 -4471.2456
[80,] -17844.3899 -16594.6620
[81,] -41669.0249 -17844.3899
[82,] -37617.2950 -41669.0249
[83,] -35739.0421 -37617.2950
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1035.1612 7394.6137
2 -8009.6633 1035.1612
3 -1979.5471 -8009.6633
4 -18938.4272 -1979.5471
5 -7691.7141 -18938.4272
6 -19109.2348 -7691.7141
7 -18728.3691 -19109.2348
8 -20245.8869 -18728.3691
9 -30697.6910 -20245.8869
10 -32842.5801 -30697.6910
11 -26146.9328 -32842.5801
12 -27422.2441 -26146.9328
13 -25850.7533 -27422.2441
14 -27860.9020 -25850.7533
15 -19203.5887 -27860.9020
16 -16505.4842 -19203.5887
17 -12009.2351 -16505.4842
18 -26725.4035 -12009.2351
19 -17006.1162 -26725.4035
20 1880.7463 -17006.1162
21 -2800.7622 1880.7463
22 -4812.0458 -2800.7622
23 3480.0805 -4812.0458
24 -9694.7793 3480.0805
25 -2933.1339 -9694.7793
26 1385.6747 -2933.1339
27 -9139.3629 1385.6747
28 -3481.2586 -9139.3629
29 716.6947 -3481.2586
30 3969.2019 716.6947
31 3163.1233 3969.2019
32 -958.2526 3163.1233
33 -1149.6629 -958.2526
34 11736.1515 -1149.6629
35 5387.5041 11736.1515
36 6319.3194 5387.5041
37 4702.7963 6319.3194
38 12975.4215 4702.7963
39 2438.8772 12975.4215
40 5660.3760 2438.8772
41 23158.7935 5660.3760
42 18958.5126 23158.7935
43 25394.8984 18958.5126
44 20335.0435 25394.8984
45 18482.8585 20335.0435
46 27981.5599 18482.8585
47 13604.4767 27981.5599
48 14871.1652 13604.4767
49 18523.0644 14871.1652
50 11197.1696 18523.0644
51 7281.2723 11197.1696
52 12889.6582 7281.2723
53 -425.3869 12889.6582
54 9382.9647 -425.3869
55 19355.3086 9382.9647
56 15927.4958 19355.3086
57 31669.7609 15927.4958
58 27284.2239 31669.7609
59 23527.8581 27284.2239
60 22745.2646 23527.8581
61 3361.1372 22745.2646
62 17477.8465 3361.1372
63 13666.4704 17477.8465
64 21044.7997 13666.4704
65 -7247.8922 21044.7997
66 17995.2046 -7247.8922
67 4415.8172 17995.2046
68 905.2438 4415.8172
69 26164.5216 905.2438
70 8269.9856 26164.5216
71 15886.0554 8269.9856
72 -14213.3394 15886.0554
73 1161.7281 -14213.3394
74 -7165.5471 1161.7281
75 6935.8788 -7165.5471
76 -669.6639 6935.8788
77 3498.7402 -669.6639
78 -4471.2456 3498.7402
79 -16594.6620 -4471.2456
80 -17844.3899 -16594.6620
81 -41669.0249 -17844.3899
82 -37617.2950 -41669.0249
83 -35739.0421 -37617.2950
> 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/7f8wm1229523647.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/8aw6x1229523647.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/9c7y01229523647.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/107j221229523647.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/11z3kl1229523647.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/1238ld1229523647.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/13g57n1229523648.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/14enqn1229523648.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/15y3rr1229523648.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/16280d1229523648.tab")
+ }
>
> system("convert tmp/1fc841229523647.ps tmp/1fc841229523647.png")
> system("convert tmp/2da4v1229523647.ps tmp/2da4v1229523647.png")
> system("convert tmp/33e3h1229523647.ps tmp/33e3h1229523647.png")
> system("convert tmp/418c31229523647.ps tmp/418c31229523647.png")
> system("convert tmp/5bp3j1229523647.ps tmp/5bp3j1229523647.png")
> system("convert tmp/65r8z1229523647.ps tmp/65r8z1229523647.png")
> system("convert tmp/7f8wm1229523647.ps tmp/7f8wm1229523647.png")
> system("convert tmp/8aw6x1229523647.ps tmp/8aw6x1229523647.png")
> system("convert tmp/9c7y01229523647.ps tmp/9c7y01229523647.png")
> system("convert tmp/107j221229523647.ps tmp/107j221229523647.png")
>
>
> proc.time()
user system elapsed
5.559 2.765 5.938