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(112.3
+ ,112.9
+ ,88.7
+ ,105.1
+ ,117.3
+ ,130.5
+ ,94.6
+ ,114.9
+ ,111.1
+ ,137.9
+ ,98.7
+ ,106.4
+ ,102.2
+ ,115
+ ,84.2
+ ,104.5
+ ,104.3
+ ,116.8
+ ,87.7
+ ,121.6
+ ,122.9
+ ,140.9
+ ,103.3
+ ,141.4
+ ,107.6
+ ,120.7
+ ,88.2
+ ,99
+ ,121.3
+ ,134.2
+ ,93.4
+ ,126.7
+ ,131.5
+ ,147.3
+ ,106.3
+ ,134.1
+ ,89
+ ,112.4
+ ,73.1
+ ,81.3
+ ,104.4
+ ,107.1
+ ,78.6
+ ,88.6
+ ,128.9
+ ,128.4
+ ,101.6
+ ,132.7
+ ,135.9
+ ,137.7
+ ,101.4
+ ,132.9
+ ,133.3
+ ,135
+ ,98.5
+ ,134.4
+ ,121.3
+ ,151
+ ,99
+ ,103.7
+ ,120.5
+ ,137.4
+ ,89.5
+ ,119.7
+ ,120.4
+ ,132.4
+ ,83.5
+ ,115
+ ,137.9
+ ,161.3
+ ,97.4
+ ,132.9
+ ,126.1
+ ,139.8
+ ,87.8
+ ,108.5
+ ,133.2
+ ,146
+ ,90.4
+ ,113.9
+ ,151.1
+ ,166.5
+ ,101.6
+ ,142
+ ,105
+ ,143.3
+ ,80
+ ,97.7
+ ,119
+ ,121
+ ,81.7
+ ,92.2
+ ,140.4
+ ,152.6
+ ,96.4
+ ,128.8
+ ,156.6
+ ,154.4
+ ,110.2
+ ,134.9
+ ,137.1
+ ,154.6
+ ,101.1
+ ,128.2
+ ,122.7
+ ,158
+ ,89.3
+ ,114.8
+ ,125.8
+ ,142.6
+ ,90
+ ,117.9
+ ,139.3
+ ,153.4
+ ,95.4
+ ,119.1
+ ,134.9
+ ,163.4
+ ,100.3
+ ,120.7
+ ,149.2
+ ,167.3
+ ,99.5
+ ,129.1
+ ,132.3
+ ,154.8
+ ,93.9
+ ,117.6
+ ,149
+ ,165.7
+ ,100.6
+ ,129.2
+ ,117.2
+ ,144.7
+ ,84.7
+ ,100
+ ,119.6
+ ,120.9
+ ,81.6
+ ,87
+ ,152
+ ,152.8
+ ,109
+ ,128
+ ,149.4
+ ,160.2
+ ,99
+ ,127.7
+ ,127.3
+ ,128.3
+ ,81.1
+ ,93.4
+ ,114.1
+ ,150.5
+ ,81.8
+ ,84.1
+ ,102.1
+ ,117
+ ,66.5
+ ,71.7
+ ,107.7
+ ,116
+ ,66.4
+ ,83.2
+ ,104.4
+ ,133.3
+ ,86.3
+ ,89.1
+ ,102.1
+ ,116.4
+ ,73.6
+ ,79.6
+ ,96
+ ,104
+ ,71.5
+ ,62.8
+ ,109.3
+ ,126.6
+ ,87.2
+ ,95.1
+ ,90
+ ,92.9
+ ,65.3
+ ,63.6
+ ,83.9
+ ,83.6
+ ,69.7
+ ,61.4
+ ,112
+ ,112.8
+ ,95.5
+ ,98.2
+ ,114.3
+ ,113.2
+ ,86.3
+ ,95.3
+ ,103.6
+ ,118.5
+ ,81
+ ,81.5
+ ,91.7
+ ,125.5
+ ,88.7
+ ,85.5
+ ,80.8
+ ,91.3
+ ,71.9
+ ,71.1
+ ,87.2
+ ,105.4
+ ,78.6
+ ,78.1
+ ,109.2
+ ,121.3
+ ,96
+ ,103
+ ,102.7
+ ,106.9
+ ,81.1
+ ,86
+ ,95.1
+ ,109.4
+ ,77.5
+ ,86.2
+ ,117.5
+ ,132.6
+ ,97.3
+ ,105.7
+ ,85.1
+ ,96.8
+ ,78.6
+ ,57.2
+ ,92.1
+ ,100.3
+ ,79
+ ,73.7
+ ,113.5
+ ,119.2
+ ,93.4
+ ,120.5)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('X1'
+ ,'X2'
+ ,'X3'
+ ,'X4')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('X1','X2','X3','X4'),1:60))
> 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 = '4'
> #'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
X4 X1 X2 X3
1 105.1 112.3 112.9 88.7
2 114.9 117.3 130.5 94.6
3 106.4 111.1 137.9 98.7
4 104.5 102.2 115.0 84.2
5 121.6 104.3 116.8 87.7
6 141.4 122.9 140.9 103.3
7 99.0 107.6 120.7 88.2
8 126.7 121.3 134.2 93.4
9 134.1 131.5 147.3 106.3
10 81.3 89.0 112.4 73.1
11 88.6 104.4 107.1 78.6
12 132.7 128.9 128.4 101.6
13 132.9 135.9 137.7 101.4
14 134.4 133.3 135.0 98.5
15 103.7 121.3 151.0 99.0
16 119.7 120.5 137.4 89.5
17 115.0 120.4 132.4 83.5
18 132.9 137.9 161.3 97.4
19 108.5 126.1 139.8 87.8
20 113.9 133.2 146.0 90.4
21 142.0 151.1 166.5 101.6
22 97.7 105.0 143.3 80.0
23 92.2 119.0 121.0 81.7
24 128.8 140.4 152.6 96.4
25 134.9 156.6 154.4 110.2
26 128.2 137.1 154.6 101.1
27 114.8 122.7 158.0 89.3
28 117.9 125.8 142.6 90.0
29 119.1 139.3 153.4 95.4
30 120.7 134.9 163.4 100.3
31 129.1 149.2 167.3 99.5
32 117.6 132.3 154.8 93.9
33 129.2 149.0 165.7 100.6
34 100.0 117.2 144.7 84.7
35 87.0 119.6 120.9 81.6
36 128.0 152.0 152.8 109.0
37 127.7 149.4 160.2 99.0
38 93.4 127.3 128.3 81.1
39 84.1 114.1 150.5 81.8
40 71.7 102.1 117.0 66.5
41 83.2 107.7 116.0 66.4
42 89.1 104.4 133.3 86.3
43 79.6 102.1 116.4 73.6
44 62.8 96.0 104.0 71.5
45 95.1 109.3 126.6 87.2
46 63.6 90.0 92.9 65.3
47 61.4 83.9 83.6 69.7
48 98.2 112.0 112.8 95.5
49 95.3 114.3 113.2 86.3
50 81.5 103.6 118.5 81.0
51 85.5 91.7 125.5 88.7
52 71.1 80.8 91.3 71.9
53 78.1 87.2 105.4 78.6
54 103.0 109.2 121.3 96.0
55 86.0 102.7 106.9 81.1
56 86.2 95.1 109.4 77.5
57 105.7 117.5 132.6 97.3
58 57.2 85.1 96.8 78.6
59 73.7 92.1 100.3 79.0
60 120.5 113.5 119.2 93.4
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X1 X2 X3
-53.4835 0.4501 0.0638 1.0923
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19.6567 -6.5793 -0.6915 5.7124 24.8845
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -53.4834 9.6830 -5.523 8.93e-07 ***
X1 0.4501 0.1475 3.052 0.00347 **
X2 0.0638 0.1190 0.536 0.59398
X3 1.0923 0.1820 6.002 1.51e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.17 on 56 degrees of freedom
Multiple R-squared: 0.844, Adjusted R-squared: 0.8356
F-statistic: 101 on 3 and 56 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.9652277 0.06954459 0.03477229
[2,] 0.9604684 0.07906327 0.03953163
[3,] 0.9334994 0.13300116 0.06650058
[4,] 0.8957606 0.20847886 0.10423943
[5,] 0.8802024 0.23959523 0.11979761
[6,] 0.8481826 0.30363490 0.15181745
[7,] 0.8026264 0.39474726 0.19737363
[8,] 0.8054544 0.38909117 0.19454559
[9,] 0.8987611 0.20247779 0.10123889
[10,] 0.9325278 0.13494439 0.06747220
[11,] 0.9521794 0.09564120 0.04782060
[12,] 0.9436928 0.11261447 0.05630723
[13,] 0.9356798 0.12864033 0.06432016
[14,] 0.9210855 0.15782904 0.07891452
[15,] 0.9107720 0.17845602 0.08922801
[16,] 0.9013584 0.19728314 0.09864157
[17,] 0.9246749 0.15065026 0.07532513
[18,] 0.9115216 0.17695678 0.08847839
[19,] 0.9388819 0.12223624 0.06111812
[20,] 0.9186082 0.16278356 0.08139178
[21,] 0.9091460 0.18170796 0.09085398
[22,] 0.9220515 0.15589697 0.07794848
[23,] 0.8991706 0.20165883 0.10082941
[24,] 0.8740286 0.25194284 0.12597142
[25,] 0.8322193 0.33556135 0.16778068
[26,] 0.7947203 0.41055936 0.20527968
[27,] 0.7420565 0.51588705 0.25794352
[28,] 0.6983285 0.60334293 0.30167147
[29,] 0.7573961 0.48520782 0.24260391
[30,] 0.8105716 0.37885672 0.18942836
[31,] 0.7491376 0.50172471 0.25086235
[32,] 0.7425949 0.51481011 0.25740506
[33,] 0.7986119 0.40277630 0.20138815
[34,] 0.7409963 0.51800740 0.25900370
[35,] 0.7073150 0.58537002 0.29268501
[36,] 0.6777827 0.64443469 0.32221735
[37,] 0.5999514 0.80009714 0.40004857
[38,] 0.6871405 0.62571893 0.31285946
[39,] 0.6153650 0.76927000 0.38463500
[40,] 0.5251036 0.94979281 0.47489640
[41,] 0.4453358 0.89067162 0.55466419
[42,] 0.4115084 0.82301682 0.58849159
[43,] 0.3424913 0.68498250 0.65750875
[44,] 0.3636422 0.72728437 0.63635782
[45,] 0.2747888 0.54957752 0.72521124
[46,] 0.2570652 0.51413037 0.74293482
[47,] 0.2731431 0.54628620 0.72685690
> postscript(file="/var/www/html/rcomp/tmp/1j9k01292933678.ps",horizontal=F,onefile=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/2j9k01292933678.ps",horizontal=F,onefile=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/3j9k01292933678.ps",horizontal=F,onefile=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/4c01l1292933678.ps",horizontal=F,onefile=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/5c01l1292933678.ps",horizontal=F,onefile=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 = 60
Frequency = 1
1 2 3 4 5
3.939932487 3.921526285 -6.738441445 12.667768515 24.884455175
6 7 8 9 10
17.733872451 0.004029655 14.995755459 2.877407530 7.700338840
11 12 13 14 15
2.398660014 8.987642928 5.661846573 11.672241232 -15.193270236
16 17 18 19 20
12.411821273 14.629910362 7.625191094 0.394963907 -0.636612754
21 22 23 24 25
5.863895785 7.389601584 -4.846330113 4.047310840 -12.333922582
26 27 28 29 30
-0.328889348 5.425649615 7.348196830 -4.116249923 -6.526210539
31 32 33 34 35
-3.937986386 -0.916172913 -4.847451780 -1.025283268 -10.200791584
36 37 38 39 40
-15.750441180 -4.428832447 -7.192746462 -12.732158589 -0.880359301
41 42 43 44 45
8.271968245 -7.184042510 -0.697710740 -11.666847993 -3.945286501
46 47 48 49 50
-0.685337615 -4.352495187 -10.246580011 -4.157839809 -7.690225994
51 52 53 54 55
-7.191382832 3.848453797 -0.250690537 -5.274729505 -2.154220433
56 57 58 59 60
5.239668250 -8.451812885 -19.656710211 -6.967849426 13.263804315
> postscript(file="/var/www/html/rcomp/tmp/6c01l1292933678.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 3.939932487 NA
1 3.921526285 3.939932487
2 -6.738441445 3.921526285
3 12.667768515 -6.738441445
4 24.884455175 12.667768515
5 17.733872451 24.884455175
6 0.004029655 17.733872451
7 14.995755459 0.004029655
8 2.877407530 14.995755459
9 7.700338840 2.877407530
10 2.398660014 7.700338840
11 8.987642928 2.398660014
12 5.661846573 8.987642928
13 11.672241232 5.661846573
14 -15.193270236 11.672241232
15 12.411821273 -15.193270236
16 14.629910362 12.411821273
17 7.625191094 14.629910362
18 0.394963907 7.625191094
19 -0.636612754 0.394963907
20 5.863895785 -0.636612754
21 7.389601584 5.863895785
22 -4.846330113 7.389601584
23 4.047310840 -4.846330113
24 -12.333922582 4.047310840
25 -0.328889348 -12.333922582
26 5.425649615 -0.328889348
27 7.348196830 5.425649615
28 -4.116249923 7.348196830
29 -6.526210539 -4.116249923
30 -3.937986386 -6.526210539
31 -0.916172913 -3.937986386
32 -4.847451780 -0.916172913
33 -1.025283268 -4.847451780
34 -10.200791584 -1.025283268
35 -15.750441180 -10.200791584
36 -4.428832447 -15.750441180
37 -7.192746462 -4.428832447
38 -12.732158589 -7.192746462
39 -0.880359301 -12.732158589
40 8.271968245 -0.880359301
41 -7.184042510 8.271968245
42 -0.697710740 -7.184042510
43 -11.666847993 -0.697710740
44 -3.945286501 -11.666847993
45 -0.685337615 -3.945286501
46 -4.352495187 -0.685337615
47 -10.246580011 -4.352495187
48 -4.157839809 -10.246580011
49 -7.690225994 -4.157839809
50 -7.191382832 -7.690225994
51 3.848453797 -7.191382832
52 -0.250690537 3.848453797
53 -5.274729505 -0.250690537
54 -2.154220433 -5.274729505
55 5.239668250 -2.154220433
56 -8.451812885 5.239668250
57 -19.656710211 -8.451812885
58 -6.967849426 -19.656710211
59 13.263804315 -6.967849426
60 NA 13.263804315
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.921526285 3.939932487
[2,] -6.738441445 3.921526285
[3,] 12.667768515 -6.738441445
[4,] 24.884455175 12.667768515
[5,] 17.733872451 24.884455175
[6,] 0.004029655 17.733872451
[7,] 14.995755459 0.004029655
[8,] 2.877407530 14.995755459
[9,] 7.700338840 2.877407530
[10,] 2.398660014 7.700338840
[11,] 8.987642928 2.398660014
[12,] 5.661846573 8.987642928
[13,] 11.672241232 5.661846573
[14,] -15.193270236 11.672241232
[15,] 12.411821273 -15.193270236
[16,] 14.629910362 12.411821273
[17,] 7.625191094 14.629910362
[18,] 0.394963907 7.625191094
[19,] -0.636612754 0.394963907
[20,] 5.863895785 -0.636612754
[21,] 7.389601584 5.863895785
[22,] -4.846330113 7.389601584
[23,] 4.047310840 -4.846330113
[24,] -12.333922582 4.047310840
[25,] -0.328889348 -12.333922582
[26,] 5.425649615 -0.328889348
[27,] 7.348196830 5.425649615
[28,] -4.116249923 7.348196830
[29,] -6.526210539 -4.116249923
[30,] -3.937986386 -6.526210539
[31,] -0.916172913 -3.937986386
[32,] -4.847451780 -0.916172913
[33,] -1.025283268 -4.847451780
[34,] -10.200791584 -1.025283268
[35,] -15.750441180 -10.200791584
[36,] -4.428832447 -15.750441180
[37,] -7.192746462 -4.428832447
[38,] -12.732158589 -7.192746462
[39,] -0.880359301 -12.732158589
[40,] 8.271968245 -0.880359301
[41,] -7.184042510 8.271968245
[42,] -0.697710740 -7.184042510
[43,] -11.666847993 -0.697710740
[44,] -3.945286501 -11.666847993
[45,] -0.685337615 -3.945286501
[46,] -4.352495187 -0.685337615
[47,] -10.246580011 -4.352495187
[48,] -4.157839809 -10.246580011
[49,] -7.690225994 -4.157839809
[50,] -7.191382832 -7.690225994
[51,] 3.848453797 -7.191382832
[52,] -0.250690537 3.848453797
[53,] -5.274729505 -0.250690537
[54,] -2.154220433 -5.274729505
[55,] 5.239668250 -2.154220433
[56,] -8.451812885 5.239668250
[57,] -19.656710211 -8.451812885
[58,] -6.967849426 -19.656710211
[59,] 13.263804315 -6.967849426
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.921526285 3.939932487
2 -6.738441445 3.921526285
3 12.667768515 -6.738441445
4 24.884455175 12.667768515
5 17.733872451 24.884455175
6 0.004029655 17.733872451
7 14.995755459 0.004029655
8 2.877407530 14.995755459
9 7.700338840 2.877407530
10 2.398660014 7.700338840
11 8.987642928 2.398660014
12 5.661846573 8.987642928
13 11.672241232 5.661846573
14 -15.193270236 11.672241232
15 12.411821273 -15.193270236
16 14.629910362 12.411821273
17 7.625191094 14.629910362
18 0.394963907 7.625191094
19 -0.636612754 0.394963907
20 5.863895785 -0.636612754
21 7.389601584 5.863895785
22 -4.846330113 7.389601584
23 4.047310840 -4.846330113
24 -12.333922582 4.047310840
25 -0.328889348 -12.333922582
26 5.425649615 -0.328889348
27 7.348196830 5.425649615
28 -4.116249923 7.348196830
29 -6.526210539 -4.116249923
30 -3.937986386 -6.526210539
31 -0.916172913 -3.937986386
32 -4.847451780 -0.916172913
33 -1.025283268 -4.847451780
34 -10.200791584 -1.025283268
35 -15.750441180 -10.200791584
36 -4.428832447 -15.750441180
37 -7.192746462 -4.428832447
38 -12.732158589 -7.192746462
39 -0.880359301 -12.732158589
40 8.271968245 -0.880359301
41 -7.184042510 8.271968245
42 -0.697710740 -7.184042510
43 -11.666847993 -0.697710740
44 -3.945286501 -11.666847993
45 -0.685337615 -3.945286501
46 -4.352495187 -0.685337615
47 -10.246580011 -4.352495187
48 -4.157839809 -10.246580011
49 -7.690225994 -4.157839809
50 -7.191382832 -7.690225994
51 3.848453797 -7.191382832
52 -0.250690537 3.848453797
53 -5.274729505 -0.250690537
54 -2.154220433 -5.274729505
55 5.239668250 -2.154220433
56 -8.451812885 5.239668250
57 -19.656710211 -8.451812885
58 -6.967849426 -19.656710211
59 13.263804315 -6.967849426
> 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/749i61292933678.ps",horizontal=F,onefile=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/8f1ir1292933678.ps",horizontal=F,onefile=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/9f1ir1292933678.ps",horizontal=F,onefile=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/10f1ir1292933678.ps",horizontal=F,onefile=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/111jyx1292933678.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/1242fk1292933678.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/13lde91292933679.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/14wmwu1292933679.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/15hmu01292933679.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/16vwsr1292933679.tab")
+ }
>
> try(system("convert tmp/1j9k01292933678.ps tmp/1j9k01292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/2j9k01292933678.ps tmp/2j9k01292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/3j9k01292933678.ps tmp/3j9k01292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/4c01l1292933678.ps tmp/4c01l1292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/5c01l1292933678.ps tmp/5c01l1292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/6c01l1292933678.ps tmp/6c01l1292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/749i61292933678.ps tmp/749i61292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/8f1ir1292933678.ps tmp/8f1ir1292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/9f1ir1292933678.ps tmp/9f1ir1292933678.png",intern=TRUE))
character(0)
> try(system("convert tmp/10f1ir1292933678.ps tmp/10f1ir1292933678.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.502 1.644 5.999