R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(106370
+ ,100.3
+ ,109375
+ ,101.9
+ ,116476
+ ,102.1
+ ,123297
+ ,103.2
+ ,114813
+ ,103.7
+ ,117925
+ ,106.2
+ ,126466
+ ,107.7
+ ,131235
+ ,109.9
+ ,120546
+ ,111.7
+ ,123791
+ ,114.9
+ ,129813
+ ,116
+ ,133463
+ ,118.3
+ ,122987
+ ,120.4
+ ,125418
+ ,126
+ ,130199
+ ,128.1
+ ,133016
+ ,130.1
+ ,121454
+ ,130.8
+ ,122044
+ ,133.6
+ ,128313
+ ,134.2
+ ,131556
+ ,135.5
+ ,120027
+ ,136.2
+ ,123001
+ ,139.1
+ ,130111
+ ,139
+ ,132524
+ ,139.6
+ ,123742
+ ,138.7
+ ,124931
+ ,140.9
+ ,133646
+ ,141.3
+ ,136557
+ ,141.8
+ ,127509
+ ,142
+ ,128945
+ ,144.5
+ ,137191
+ ,144.6
+ ,139716
+ ,145.5
+ ,129083
+ ,146.8
+ ,131604
+ ,149.5
+ ,139413
+ ,149.9
+ ,143125
+ ,150.1
+ ,133948
+ ,150.9
+ ,137116
+ ,152.8
+ ,144864
+ ,153.1
+ ,149277
+ ,154
+ ,138796
+ ,154.9
+ ,143258
+ ,156.9
+ ,150034
+ ,158.4
+ ,154708
+ ,159.7
+ ,144888
+ ,160.2
+ ,148762
+ ,163.2
+ ,156500
+ ,163.7
+ ,161088
+ ,164.4
+ ,152772
+ ,163.7
+ ,158011
+ ,165.5
+ ,163318
+ ,165.6
+ ,169969
+ ,166.8
+ ,162269
+ ,167.5
+ ,165765
+ ,170.6
+ ,170600
+ ,170.9
+ ,174681
+ ,172
+ ,166364
+ ,171.8
+ ,170240
+ ,173.9
+ ,176150
+ ,174
+ ,182056
+ ,173.8
+ ,172218
+ ,173.9
+ ,177856
+ ,176
+ ,182253
+ ,176.6
+ ,188090
+ ,178.2
+ ,176863
+ ,179.2
+ ,183273
+ ,181.3
+ ,187969
+ ,181.8
+ ,194650
+ ,182.9
+ ,183036
+ ,183.8
+ ,189516
+ ,186.3
+ ,193805
+ ,187.4
+ ,200499
+ ,189.2
+ ,188142
+ ,189.7
+ ,193732
+ ,191.9
+ ,197126
+ ,192.6
+ ,205140
+ ,193.7
+ ,191751
+ ,194.2
+ ,196700
+ ,197.6
+ ,199784
+ ,199.3
+ ,207360
+ ,201.4
+ ,196101
+ ,203
+ ,200824
+ ,206.3
+ ,205743
+ ,207.1
+ ,212489
+ ,209.8
+ ,200810
+ ,211.1
+ ,203683
+ ,215.3
+ ,207286
+ ,217.4
+ ,210910
+ ,215.5
+ ,194915
+ ,210.9
+ ,217920
+ ,212.6)
+ ,dim=c(2
+ ,90)
+ ,dimnames=list(c('HFCE'
+ ,'RPI')
+ ,1:90))
> y <- array(NA,dim=c(2,90),dimnames=list(c('HFCE','RPI'),1:90))
> 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 Quarterly 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
HFCE RPI Q1 Q2 Q3 t
1 106370 100.3 1 0 0 1
2 109375 101.9 0 1 0 2
3 116476 102.1 0 0 1 3
4 123297 103.2 0 0 0 4
5 114813 103.7 1 0 0 5
6 117925 106.2 0 1 0 6
7 126466 107.7 0 0 1 7
8 131235 109.9 0 0 0 8
9 120546 111.7 1 0 0 9
10 123791 114.9 0 1 0 10
11 129813 116.0 0 0 1 11
12 133463 118.3 0 0 0 12
13 122987 120.4 1 0 0 13
14 125418 126.0 0 1 0 14
15 130199 128.1 0 0 1 15
16 133016 130.1 0 0 0 16
17 121454 130.8 1 0 0 17
18 122044 133.6 0 1 0 18
19 128313 134.2 0 0 1 19
20 131556 135.5 0 0 0 20
21 120027 136.2 1 0 0 21
22 123001 139.1 0 1 0 22
23 130111 139.0 0 0 1 23
24 132524 139.6 0 0 0 24
25 123742 138.7 1 0 0 25
26 124931 140.9 0 1 0 26
27 133646 141.3 0 0 1 27
28 136557 141.8 0 0 0 28
29 127509 142.0 1 0 0 29
30 128945 144.5 0 1 0 30
31 137191 144.6 0 0 1 31
32 139716 145.5 0 0 0 32
33 129083 146.8 1 0 0 33
34 131604 149.5 0 1 0 34
35 139413 149.9 0 0 1 35
36 143125 150.1 0 0 0 36
37 133948 150.9 1 0 0 37
38 137116 152.8 0 1 0 38
39 144864 153.1 0 0 1 39
40 149277 154.0 0 0 0 40
41 138796 154.9 1 0 0 41
42 143258 156.9 0 1 0 42
43 150034 158.4 0 0 1 43
44 154708 159.7 0 0 0 44
45 144888 160.2 1 0 0 45
46 148762 163.2 0 1 0 46
47 156500 163.7 0 0 1 47
48 161088 164.4 0 0 0 48
49 152772 163.7 1 0 0 49
50 158011 165.5 0 1 0 50
51 163318 165.6 0 0 1 51
52 169969 166.8 0 0 0 52
53 162269 167.5 1 0 0 53
54 165765 170.6 0 1 0 54
55 170600 170.9 0 0 1 55
56 174681 172.0 0 0 0 56
57 166364 171.8 1 0 0 57
58 170240 173.9 0 1 0 58
59 176150 174.0 0 0 1 59
60 182056 173.8 0 0 0 60
61 172218 173.9 1 0 0 61
62 177856 176.0 0 1 0 62
63 182253 176.6 0 0 1 63
64 188090 178.2 0 0 0 64
65 176863 179.2 1 0 0 65
66 183273 181.3 0 1 0 66
67 187969 181.8 0 0 1 67
68 194650 182.9 0 0 0 68
69 183036 183.8 1 0 0 69
70 189516 186.3 0 1 0 70
71 193805 187.4 0 0 1 71
72 200499 189.2 0 0 0 72
73 188142 189.7 1 0 0 73
74 193732 191.9 0 1 0 74
75 197126 192.6 0 0 1 75
76 205140 193.7 0 0 0 76
77 191751 194.2 1 0 0 77
78 196700 197.6 0 1 0 78
79 199784 199.3 0 0 1 79
80 207360 201.4 0 0 0 80
81 196101 203.0 1 0 0 81
82 200824 206.3 0 1 0 82
83 205743 207.1 0 0 1 83
84 212489 209.8 0 0 0 84
85 200810 211.1 1 0 0 85
86 203683 215.3 0 1 0 86
87 207286 217.4 0 0 1 87
88 210910 215.5 0 0 0 88
89 194915 210.9 1 0 0 89
90 217920 212.6 0 1 0 90
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) RPI Q1 Q2 Q3 t
207810.6 -916.4 -12387.1 -7662.9 -3693.2 2256.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8746.32 -3560.34 83.04 3541.61 9628.35
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 207810.6 15141.0 13.725 < 2e-16 ***
RPI -916.4 142.8 -6.419 7.78e-09 ***
Q1 -12387.1 1503.8 -8.237 2.03e-12 ***
Q2 -7662.9 1498.3 -5.115 1.95e-06 ***
Q3 -3693.2 1513.5 -2.440 0.0168 *
t 2256.4 172.3 13.096 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5019 on 84 degrees of freedom
Multiple R-squared: 0.9755, Adjusted R-squared: 0.9741
F-statistic: 670.1 on 5 and 84 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.0133604221 0.0267208441 0.986639578
[2,] 0.0025848901 0.0051697801 0.997415110
[3,] 0.0010774086 0.0021548171 0.998922591
[4,] 0.0010759049 0.0021518099 0.998924095
[5,] 0.0007535473 0.0015070946 0.999246453
[6,] 0.0004069355 0.0008138709 0.999593065
[7,] 0.0002007400 0.0004014800 0.999799260
[8,] 0.0002828884 0.0005657769 0.999717112
[9,] 0.0037691558 0.0075383117 0.996230844
[10,] 0.0468723966 0.0937447932 0.953127603
[11,] 0.1887962429 0.3775924857 0.811203757
[12,] 0.4869445262 0.9738890523 0.513055474
[13,] 0.7823894266 0.4352211467 0.217610573
[14,] 0.8773258635 0.2453482730 0.122674136
[15,] 0.9440848908 0.1118302183 0.055915109
[16,] 0.9758492673 0.0483014654 0.024150733
[17,] 0.9851139452 0.0297721097 0.014886055
[18,] 0.9837812817 0.0324374366 0.016218718
[19,] 0.9870774576 0.0258450848 0.012922542
[20,] 0.9845788775 0.0308422450 0.015421122
[21,] 0.9867557186 0.0264885627 0.013244281
[22,] 0.9801699200 0.0396601600 0.019830080
[23,] 0.9785723325 0.0428553351 0.021427668
[24,] 0.9687580235 0.0624839530 0.031241976
[25,] 0.9578668985 0.0842662029 0.042133101
[26,] 0.9432829811 0.1134340379 0.056717019
[27,] 0.9245360428 0.1509279145 0.075463957
[28,] 0.9041995992 0.1916008017 0.095800401
[29,] 0.8837851648 0.2324296705 0.116214835
[30,] 0.8854720918 0.2290558163 0.114527908
[31,] 0.8684885782 0.2630228436 0.131511422
[32,] 0.8656585681 0.2686828637 0.134341432
[33,] 0.8523286129 0.2953427743 0.147671387
[34,] 0.8907703363 0.2184593273 0.109229664
[35,] 0.8861824442 0.2276351115 0.113817556
[36,] 0.8987232488 0.2025535023 0.101276751
[37,] 0.9028014315 0.1943971369 0.097198568
[38,] 0.9392059424 0.1215881152 0.060794058
[39,] 0.9456326655 0.1087346689 0.054367334
[40,] 0.9597038445 0.0805923111 0.040296156
[41,] 0.9676684202 0.0646631597 0.032331580
[42,] 0.9831109490 0.0337781021 0.016889051
[43,] 0.9839735528 0.0320528945 0.016026447
[44,] 0.9878614266 0.0242771468 0.012138573
[45,] 0.9935748524 0.0128502952 0.006425148
[46,] 0.9955516579 0.0088966842 0.004448342
[47,] 0.9948341164 0.0103317673 0.005165884
[48,] 0.9946205994 0.0107588012 0.005379401
[49,] 0.9938517705 0.0122964590 0.006148230
[50,] 0.9942138701 0.0115722598 0.005786130
[51,] 0.9919827537 0.0160344926 0.008017246
[52,] 0.9906813796 0.0186372408 0.009318620
[53,] 0.9861916454 0.0276167092 0.013808355
[54,] 0.9848109364 0.0303781271 0.015189064
[55,] 0.9784663082 0.0430673837 0.021533692
[56,] 0.9742013676 0.0515972649 0.025798632
[57,] 0.9627290396 0.0745419208 0.037270960
[58,] 0.9606348877 0.0787302247 0.039365112
[59,] 0.9447543777 0.1104912446 0.055245622
[60,] 0.9284640438 0.1430719124 0.071535956
[61,] 0.8991121376 0.2017757248 0.100887862
[62,] 0.8800620893 0.2398758213 0.119937911
[63,] 0.8373408676 0.3253182647 0.162659132
[64,] 0.7987524310 0.4024951380 0.201247569
[65,] 0.7441095260 0.5117809479 0.255890474
[66,] 0.6901728537 0.6196542926 0.309827146
[67,] 0.5998476997 0.8003046007 0.400152300
[68,] 0.5207585625 0.9584828750 0.479241437
[69,] 0.4224703663 0.8449407326 0.577529634
[70,] 0.3699464444 0.7398928889 0.630053556
[71,] 0.2908649932 0.5817299864 0.709135007
[72,] 0.2000780680 0.4001561359 0.799921932
[73,] 0.1312748641 0.2625497283 0.868725136
> postscript(file="/var/www/html/rcomp/tmp/15pf31258725244.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/2icme1258725244.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/3y2421258725244.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/44b151258725244.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/557vq1258725244.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 = 90
Frequency = 1
1 2 3 4 5 6
606.41596 -1902.92140 -844.78557 1034.73420 3139.68515 1562.11947
7 8 9 10 11 12
5251.59218 6087.16622 7178.45404 6375.37744 7179.28496 6987.50030
13 14 15 16 17 18
8566.71201 9149.02657 9628.34706 8328.63850 7538.87205 3714.23026
19 20 21 22 23 24
4306.93129 2791.73365 2034.96720 685.96671 1478.17865 -1508.50807
25 26 27 28 29 30
-984.53529 -4760.02486 -1904.60643 -4484.93445 -3218.90739 -6472.47307
31 32 33 34 35 36
-4360.97853 -6960.74137 -6271.66003 -8256.94312 -6307.52468 -8361.77660
37 38 39 40 41 42
-6674.90175 -8746.31522 -6949.53809 -7661.30092 -7186.78478 -7872.55694
43 44 45 46 47 48
-5948.08424 -6032.28188 -5263.33093 -5620.69011 -3650.63038 -4370.67581
49 50 51 52 53 54
-3197.42043 -3289.47520 -4116.98066 -2315.81960 756.41395 112.69606
55 56 57 58 59 60
-1003.52681 -1864.00704 -233.54518 -1413.67605 -1638.18151 -1864.99862
61 62 63 64 65 66
-1480.61286 -898.74373 -2178.04270 -824.31644 -1004.15901 349.71012
67 68 69 70 71 72
-722.23014 1017.28962 358.80576 2149.24008 1220.14760 3614.15645
73 74 75 76 77 78
1846.10740 2471.61783 280.96015 3353.47992 553.43087 1637.63687
79 80 81 82 83 84
53.39217 3604.32491 3942.33014 4708.89484 4134.87847 7405.65899
85 86 87 88 89 90
7048.74033 6790.07671 6091.39720 2024.67804 -8055.07720 9527.22674
> postscript(file="/var/www/html/rcomp/tmp/6u6k11258725244.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 = 90
Frequency = 1
lag(myerror, k = 1) myerror
0 606.41596 NA
1 -1902.92140 606.41596
2 -844.78557 -1902.92140
3 1034.73420 -844.78557
4 3139.68515 1034.73420
5 1562.11947 3139.68515
6 5251.59218 1562.11947
7 6087.16622 5251.59218
8 7178.45404 6087.16622
9 6375.37744 7178.45404
10 7179.28496 6375.37744
11 6987.50030 7179.28496
12 8566.71201 6987.50030
13 9149.02657 8566.71201
14 9628.34706 9149.02657
15 8328.63850 9628.34706
16 7538.87205 8328.63850
17 3714.23026 7538.87205
18 4306.93129 3714.23026
19 2791.73365 4306.93129
20 2034.96720 2791.73365
21 685.96671 2034.96720
22 1478.17865 685.96671
23 -1508.50807 1478.17865
24 -984.53529 -1508.50807
25 -4760.02486 -984.53529
26 -1904.60643 -4760.02486
27 -4484.93445 -1904.60643
28 -3218.90739 -4484.93445
29 -6472.47307 -3218.90739
30 -4360.97853 -6472.47307
31 -6960.74137 -4360.97853
32 -6271.66003 -6960.74137
33 -8256.94312 -6271.66003
34 -6307.52468 -8256.94312
35 -8361.77660 -6307.52468
36 -6674.90175 -8361.77660
37 -8746.31522 -6674.90175
38 -6949.53809 -8746.31522
39 -7661.30092 -6949.53809
40 -7186.78478 -7661.30092
41 -7872.55694 -7186.78478
42 -5948.08424 -7872.55694
43 -6032.28188 -5948.08424
44 -5263.33093 -6032.28188
45 -5620.69011 -5263.33093
46 -3650.63038 -5620.69011
47 -4370.67581 -3650.63038
48 -3197.42043 -4370.67581
49 -3289.47520 -3197.42043
50 -4116.98066 -3289.47520
51 -2315.81960 -4116.98066
52 756.41395 -2315.81960
53 112.69606 756.41395
54 -1003.52681 112.69606
55 -1864.00704 -1003.52681
56 -233.54518 -1864.00704
57 -1413.67605 -233.54518
58 -1638.18151 -1413.67605
59 -1864.99862 -1638.18151
60 -1480.61286 -1864.99862
61 -898.74373 -1480.61286
62 -2178.04270 -898.74373
63 -824.31644 -2178.04270
64 -1004.15901 -824.31644
65 349.71012 -1004.15901
66 -722.23014 349.71012
67 1017.28962 -722.23014
68 358.80576 1017.28962
69 2149.24008 358.80576
70 1220.14760 2149.24008
71 3614.15645 1220.14760
72 1846.10740 3614.15645
73 2471.61783 1846.10740
74 280.96015 2471.61783
75 3353.47992 280.96015
76 553.43087 3353.47992
77 1637.63687 553.43087
78 53.39217 1637.63687
79 3604.32491 53.39217
80 3942.33014 3604.32491
81 4708.89484 3942.33014
82 4134.87847 4708.89484
83 7405.65899 4134.87847
84 7048.74033 7405.65899
85 6790.07671 7048.74033
86 6091.39720 6790.07671
87 2024.67804 6091.39720
88 -8055.07720 2024.67804
89 9527.22674 -8055.07720
90 NA 9527.22674
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1902.92140 606.41596
[2,] -844.78557 -1902.92140
[3,] 1034.73420 -844.78557
[4,] 3139.68515 1034.73420
[5,] 1562.11947 3139.68515
[6,] 5251.59218 1562.11947
[7,] 6087.16622 5251.59218
[8,] 7178.45404 6087.16622
[9,] 6375.37744 7178.45404
[10,] 7179.28496 6375.37744
[11,] 6987.50030 7179.28496
[12,] 8566.71201 6987.50030
[13,] 9149.02657 8566.71201
[14,] 9628.34706 9149.02657
[15,] 8328.63850 9628.34706
[16,] 7538.87205 8328.63850
[17,] 3714.23026 7538.87205
[18,] 4306.93129 3714.23026
[19,] 2791.73365 4306.93129
[20,] 2034.96720 2791.73365
[21,] 685.96671 2034.96720
[22,] 1478.17865 685.96671
[23,] -1508.50807 1478.17865
[24,] -984.53529 -1508.50807
[25,] -4760.02486 -984.53529
[26,] -1904.60643 -4760.02486
[27,] -4484.93445 -1904.60643
[28,] -3218.90739 -4484.93445
[29,] -6472.47307 -3218.90739
[30,] -4360.97853 -6472.47307
[31,] -6960.74137 -4360.97853
[32,] -6271.66003 -6960.74137
[33,] -8256.94312 -6271.66003
[34,] -6307.52468 -8256.94312
[35,] -8361.77660 -6307.52468
[36,] -6674.90175 -8361.77660
[37,] -8746.31522 -6674.90175
[38,] -6949.53809 -8746.31522
[39,] -7661.30092 -6949.53809
[40,] -7186.78478 -7661.30092
[41,] -7872.55694 -7186.78478
[42,] -5948.08424 -7872.55694
[43,] -6032.28188 -5948.08424
[44,] -5263.33093 -6032.28188
[45,] -5620.69011 -5263.33093
[46,] -3650.63038 -5620.69011
[47,] -4370.67581 -3650.63038
[48,] -3197.42043 -4370.67581
[49,] -3289.47520 -3197.42043
[50,] -4116.98066 -3289.47520
[51,] -2315.81960 -4116.98066
[52,] 756.41395 -2315.81960
[53,] 112.69606 756.41395
[54,] -1003.52681 112.69606
[55,] -1864.00704 -1003.52681
[56,] -233.54518 -1864.00704
[57,] -1413.67605 -233.54518
[58,] -1638.18151 -1413.67605
[59,] -1864.99862 -1638.18151
[60,] -1480.61286 -1864.99862
[61,] -898.74373 -1480.61286
[62,] -2178.04270 -898.74373
[63,] -824.31644 -2178.04270
[64,] -1004.15901 -824.31644
[65,] 349.71012 -1004.15901
[66,] -722.23014 349.71012
[67,] 1017.28962 -722.23014
[68,] 358.80576 1017.28962
[69,] 2149.24008 358.80576
[70,] 1220.14760 2149.24008
[71,] 3614.15645 1220.14760
[72,] 1846.10740 3614.15645
[73,] 2471.61783 1846.10740
[74,] 280.96015 2471.61783
[75,] 3353.47992 280.96015
[76,] 553.43087 3353.47992
[77,] 1637.63687 553.43087
[78,] 53.39217 1637.63687
[79,] 3604.32491 53.39217
[80,] 3942.33014 3604.32491
[81,] 4708.89484 3942.33014
[82,] 4134.87847 4708.89484
[83,] 7405.65899 4134.87847
[84,] 7048.74033 7405.65899
[85,] 6790.07671 7048.74033
[86,] 6091.39720 6790.07671
[87,] 2024.67804 6091.39720
[88,] -8055.07720 2024.67804
[89,] 9527.22674 -8055.07720
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1902.92140 606.41596
2 -844.78557 -1902.92140
3 1034.73420 -844.78557
4 3139.68515 1034.73420
5 1562.11947 3139.68515
6 5251.59218 1562.11947
7 6087.16622 5251.59218
8 7178.45404 6087.16622
9 6375.37744 7178.45404
10 7179.28496 6375.37744
11 6987.50030 7179.28496
12 8566.71201 6987.50030
13 9149.02657 8566.71201
14 9628.34706 9149.02657
15 8328.63850 9628.34706
16 7538.87205 8328.63850
17 3714.23026 7538.87205
18 4306.93129 3714.23026
19 2791.73365 4306.93129
20 2034.96720 2791.73365
21 685.96671 2034.96720
22 1478.17865 685.96671
23 -1508.50807 1478.17865
24 -984.53529 -1508.50807
25 -4760.02486 -984.53529
26 -1904.60643 -4760.02486
27 -4484.93445 -1904.60643
28 -3218.90739 -4484.93445
29 -6472.47307 -3218.90739
30 -4360.97853 -6472.47307
31 -6960.74137 -4360.97853
32 -6271.66003 -6960.74137
33 -8256.94312 -6271.66003
34 -6307.52468 -8256.94312
35 -8361.77660 -6307.52468
36 -6674.90175 -8361.77660
37 -8746.31522 -6674.90175
38 -6949.53809 -8746.31522
39 -7661.30092 -6949.53809
40 -7186.78478 -7661.30092
41 -7872.55694 -7186.78478
42 -5948.08424 -7872.55694
43 -6032.28188 -5948.08424
44 -5263.33093 -6032.28188
45 -5620.69011 -5263.33093
46 -3650.63038 -5620.69011
47 -4370.67581 -3650.63038
48 -3197.42043 -4370.67581
49 -3289.47520 -3197.42043
50 -4116.98066 -3289.47520
51 -2315.81960 -4116.98066
52 756.41395 -2315.81960
53 112.69606 756.41395
54 -1003.52681 112.69606
55 -1864.00704 -1003.52681
56 -233.54518 -1864.00704
57 -1413.67605 -233.54518
58 -1638.18151 -1413.67605
59 -1864.99862 -1638.18151
60 -1480.61286 -1864.99862
61 -898.74373 -1480.61286
62 -2178.04270 -898.74373
63 -824.31644 -2178.04270
64 -1004.15901 -824.31644
65 349.71012 -1004.15901
66 -722.23014 349.71012
67 1017.28962 -722.23014
68 358.80576 1017.28962
69 2149.24008 358.80576
70 1220.14760 2149.24008
71 3614.15645 1220.14760
72 1846.10740 3614.15645
73 2471.61783 1846.10740
74 280.96015 2471.61783
75 3353.47992 280.96015
76 553.43087 3353.47992
77 1637.63687 553.43087
78 53.39217 1637.63687
79 3604.32491 53.39217
80 3942.33014 3604.32491
81 4708.89484 3942.33014
82 4134.87847 4708.89484
83 7405.65899 4134.87847
84 7048.74033 7405.65899
85 6790.07671 7048.74033
86 6091.39720 6790.07671
87 2024.67804 6091.39720
88 -8055.07720 2024.67804
89 9527.22674 -8055.07720
> 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/7xf6n1258725244.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/81aqo1258725244.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/9vtnw1258725244.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/10ebvq1258725244.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/11teqz1258725244.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/12hx641258725244.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/13aukw1258725244.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/14vgee1258725244.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/15tgh11258725244.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/16lq7m1258725244.tab")
+ }
>
> system("convert tmp/15pf31258725244.ps tmp/15pf31258725244.png")
> system("convert tmp/2icme1258725244.ps tmp/2icme1258725244.png")
> system("convert tmp/3y2421258725244.ps tmp/3y2421258725244.png")
> system("convert tmp/44b151258725244.ps tmp/44b151258725244.png")
> system("convert tmp/557vq1258725244.ps tmp/557vq1258725244.png")
> system("convert tmp/6u6k11258725244.ps tmp/6u6k11258725244.png")
> system("convert tmp/7xf6n1258725244.ps tmp/7xf6n1258725244.png")
> system("convert tmp/81aqo1258725244.ps tmp/81aqo1258725244.png")
> system("convert tmp/9vtnw1258725244.ps tmp/9vtnw1258725244.png")
> system("convert tmp/10ebvq1258725244.ps tmp/10ebvq1258725244.png")
>
>
> proc.time()
user system elapsed
2.729 1.681 3.249