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(100
+ ,100
+ ,97.82226485
+ ,99.87129987
+ ,94.04971502
+ ,99.54459954
+ ,91.12460521
+ ,99.81189981
+ ,93.13202153
+ ,100.4851005
+ ,93.88342812
+ ,101.1385011
+ ,92.55349954
+ ,101.3662014
+ ,94.43494835
+ ,101.5147015
+ ,96.25017563
+ ,101.8216018
+ ,100.4355715
+ ,102.4354024
+ ,101.5036685
+ ,102.5344025
+ ,99.39789728
+ ,102.6532027
+ ,99.68990733
+ ,102.4651025
+ ,101.6895041
+ ,102.4354024
+ ,103.6652759
+ ,102.4156024
+ ,103.0532766
+ ,102.4453024
+ ,100.9500712
+ ,102.8908029
+ ,102.345366
+ ,102.8512029
+ ,101.6472299
+ ,103.3561034
+ ,99.56809393
+ ,103.7422037
+ ,95.67727392
+ ,103.7224037
+ ,96.58494865
+ ,104.0788041
+ ,96.32604937
+ ,104.2075042
+ ,95.37109101
+ ,103.9105039
+ ,96.00056203
+ ,103.7026037
+ ,96.88367859
+ ,103.960004
+ ,94.85280372
+ ,104.0986041
+ ,92.46943974
+ ,104.1481041
+ ,93.99180173
+ ,104.7124047
+ ,93.45262168
+ ,104.7223047
+ ,92.26698759
+ ,105.1975052
+ ,90.39653498
+ ,105.0688051
+ ,90.43001228
+ ,105.0589051
+ ,91.04995327
+ ,105.5044055
+ ,89.07845784
+ ,105.3757054
+ ,89.69314509
+ ,105.4747055
+ ,87.92459054
+ ,106.029106
+ ,85.8789319
+ ,107.019107
+ ,83.20612366
+ ,107.3161073
+ ,83.85722053
+ ,107.7517078
+ ,83.01393462
+ ,108.5239085
+ ,82.84508195
+ ,109.3159093
+ ,78.68864276
+ ,109.5634096
+ ,77.56959675
+ ,110.5435105
+ ,78.53689529
+ ,111.1573112
+ ,78.55717715
+ ,111.7414117
+ ,77.4761291
+ ,111.0583111
+ ,81.58931659
+ ,111.2365112
+ ,85.02428326
+ ,111.038511
+ ,91.71290159
+ ,110.3752104
+ ,95.96293061
+ ,110.1376101
+ ,90.84689022
+ ,110.2465102
+ ,92.28788036
+ ,110.6227106
+ ,95.56511274
+ ,109.98911
+ ,93.62452884
+ ,110.2168102
+ ,92.63071726
+ ,110.1376101
+ ,89.50914211
+ ,109.9297099
+ ,87.17171779
+ ,109.8604099
+ ,86.72624975
+ ,110.1970102
+ ,85.63212844
+ ,109.9099099)
+ ,dim=c(2
+ ,60)
+ ,dimnames=list(c('wisselkoers'
+ ,'consumptieprijzen')
+ ,1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('wisselkoers','consumptieprijzen'),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 = '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
wisselkoers consumptieprijzen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 100.00000 100.0000 1 0 0 0 0 0 0 0 0 0 0 1
2 97.82226 99.8713 0 1 0 0 0 0 0 0 0 0 0 2
3 94.04972 99.5446 0 0 1 0 0 0 0 0 0 0 0 3
4 91.12461 99.8119 0 0 0 1 0 0 0 0 0 0 0 4
5 93.13202 100.4851 0 0 0 0 1 0 0 0 0 0 0 5
6 93.88343 101.1385 0 0 0 0 0 1 0 0 0 0 0 6
7 92.55350 101.3662 0 0 0 0 0 0 1 0 0 0 0 7
8 94.43495 101.5147 0 0 0 0 0 0 0 1 0 0 0 8
9 96.25018 101.8216 0 0 0 0 0 0 0 0 1 0 0 9
10 100.43557 102.4354 0 0 0 0 0 0 0 0 0 1 0 10
11 101.50367 102.5344 0 0 0 0 0 0 0 0 0 0 1 11
12 99.39790 102.6532 0 0 0 0 0 0 0 0 0 0 0 12
13 99.68991 102.4651 1 0 0 0 0 0 0 0 0 0 0 13
14 101.68950 102.4354 0 1 0 0 0 0 0 0 0 0 0 14
15 103.66528 102.4156 0 0 1 0 0 0 0 0 0 0 0 15
16 103.05328 102.4453 0 0 0 1 0 0 0 0 0 0 0 16
17 100.95007 102.8908 0 0 0 0 1 0 0 0 0 0 0 17
18 102.34537 102.8512 0 0 0 0 0 1 0 0 0 0 0 18
19 101.64723 103.3561 0 0 0 0 0 0 1 0 0 0 0 19
20 99.56809 103.7422 0 0 0 0 0 0 0 1 0 0 0 20
21 95.67727 103.7224 0 0 0 0 0 0 0 0 1 0 0 21
22 96.58495 104.0788 0 0 0 0 0 0 0 0 0 1 0 22
23 96.32605 104.2075 0 0 0 0 0 0 0 0 0 0 1 23
24 95.37109 103.9105 0 0 0 0 0 0 0 0 0 0 0 24
25 96.00056 103.7026 1 0 0 0 0 0 0 0 0 0 0 25
26 96.88368 103.9600 0 1 0 0 0 0 0 0 0 0 0 26
27 94.85280 104.0986 0 0 1 0 0 0 0 0 0 0 0 27
28 92.46944 104.1481 0 0 0 1 0 0 0 0 0 0 0 28
29 93.99180 104.7124 0 0 0 0 1 0 0 0 0 0 0 29
30 93.45262 104.7223 0 0 0 0 0 1 0 0 0 0 0 30
31 92.26699 105.1975 0 0 0 0 0 0 1 0 0 0 0 31
32 90.39653 105.0688 0 0 0 0 0 0 0 1 0 0 0 32
33 90.43001 105.0589 0 0 0 0 0 0 0 0 1 0 0 33
34 91.04995 105.5044 0 0 0 0 0 0 0 0 0 1 0 34
35 89.07846 105.3757 0 0 0 0 0 0 0 0 0 0 1 35
36 89.69315 105.4747 0 0 0 0 0 0 0 0 0 0 0 36
37 87.92459 106.0291 1 0 0 0 0 0 0 0 0 0 0 37
38 85.87893 107.0191 0 1 0 0 0 0 0 0 0 0 0 38
39 83.20612 107.3161 0 0 1 0 0 0 0 0 0 0 0 39
40 83.85722 107.7517 0 0 0 1 0 0 0 0 0 0 0 40
41 83.01393 108.5239 0 0 0 0 1 0 0 0 0 0 0 41
42 82.84508 109.3159 0 0 0 0 0 1 0 0 0 0 0 42
43 78.68864 109.5634 0 0 0 0 0 0 1 0 0 0 0 43
44 77.56960 110.5435 0 0 0 0 0 0 0 1 0 0 0 44
45 78.53690 111.1573 0 0 0 0 0 0 0 0 1 0 0 45
46 78.55718 111.7414 0 0 0 0 0 0 0 0 0 1 0 46
47 77.47613 111.0583 0 0 0 0 0 0 0 0 0 0 1 47
48 81.58932 111.2365 0 0 0 0 0 0 0 0 0 0 0 48
49 85.02428 111.0385 1 0 0 0 0 0 0 0 0 0 0 49
50 91.71290 110.3752 0 1 0 0 0 0 0 0 0 0 0 50
51 95.96293 110.1376 0 0 1 0 0 0 0 0 0 0 0 51
52 90.84689 110.2465 0 0 0 1 0 0 0 0 0 0 0 52
53 92.28788 110.6227 0 0 0 0 1 0 0 0 0 0 0 53
54 95.56511 109.9891 0 0 0 0 0 1 0 0 0 0 0 54
55 93.62453 110.2168 0 0 0 0 0 0 1 0 0 0 0 55
56 92.63072 110.1376 0 0 0 0 0 0 0 1 0 0 0 56
57 89.50914 109.9297 0 0 0 0 0 0 0 0 1 0 0 57
58 87.17172 109.8604 0 0 0 0 0 0 0 0 0 1 0 58
59 86.72625 110.1970 0 0 0 0 0 0 0 0 0 0 1 59
60 85.63213 109.9099 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) consumptieprijzen M1 M2
349.9882 -2.5248 1.2949 2.3133
M3 M4 M5 M6
1.5220 -0.3713 1.1971 2.2691
M7 M8 M9 M10
0.9906 0.5481 -0.2124 1.1754
M11 t
0.2465 0.2662
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.7365 -3.7731 -0.7642 4.8242 8.9501
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 349.9882 69.0012 5.072 6.9e-06 ***
consumptieprijzen -2.5248 0.6940 -3.638 0.000693 ***
M1 1.2949 3.3664 0.385 0.702266
M2 2.3133 3.3583 0.689 0.494399
M3 1.5220 3.3536 0.454 0.652071
M4 -0.3713 3.3500 -0.111 0.912232
M5 1.1971 3.3488 0.357 0.722367
M6 2.2691 3.3443 0.678 0.500862
M7 0.9906 3.3465 0.296 0.768551
M8 0.5481 3.3472 0.164 0.870643
M9 -0.2124 3.3423 -0.064 0.949608
M10 1.1754 3.3517 0.351 0.727414
M11 0.2465 3.3377 0.074 0.941438
t 0.2662 0.1462 1.821 0.075172 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.27 on 46 degrees of freedom
Multiple R-squared: 0.5506, Adjusted R-squared: 0.4236
F-statistic: 4.335 on 13 and 46 DF, p-value: 0.0001029
> 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.27763231 0.55526462 0.72236769
[2,] 0.23755077 0.47510155 0.76244923
[3,] 0.15690623 0.31381246 0.84309377
[4,] 0.10663123 0.21326246 0.89336877
[5,] 0.13519025 0.27038051 0.86480975
[6,] 0.17906666 0.35813332 0.82093334
[7,] 0.25168850 0.50337699 0.74831150
[8,] 0.23322277 0.46644554 0.76677723
[9,] 0.20322907 0.40645813 0.79677093
[10,] 0.19043415 0.38086829 0.80956585
[11,] 0.23861533 0.47723066 0.76138467
[12,] 0.23743408 0.47486815 0.76256592
[13,] 0.20305121 0.40610241 0.79694879
[14,] 0.14861197 0.29722393 0.85138803
[15,] 0.12579288 0.25158575 0.87420712
[16,] 0.08610373 0.17220746 0.91389627
[17,] 0.06132479 0.12264958 0.93867521
[18,] 0.06177954 0.12355907 0.93822046
[19,] 0.08429559 0.16859119 0.91570441
[20,] 0.34234303 0.68468606 0.65765697
[21,] 0.77256395 0.45487210 0.22743605
[22,] 0.93442220 0.13115560 0.06557780
[23,] 0.93674725 0.12650549 0.06325275
[24,] 0.93998890 0.12002220 0.06001110
[25,] 0.96364552 0.07270897 0.03635448
[26,] 0.92903646 0.14192708 0.07096354
[27,] 0.84929196 0.30141607 0.15070804
> postscript(file="/var/www/html/rcomp/tmp/14ygz1258735091.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/2p79v1258735091.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/3179d1258735091.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/4y3il1258735091.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/5d8dx1258735091.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 = 60
Frequency = 1
1 2 3 4 5 6 7
0.9272595 -2.8599596 -6.9322791 -7.5553488 -5.6828397 -4.6198542 -4.3625855
8 9 10 11 12 13 14
-1.9298826 1.1545405 5.2356523 7.2164400 5.3909885 3.6470343 4.2870989
15 16 17 18 19 20 21
6.7379524 7.8281083 5.0151061 4.9723158 6.5612424 5.6332454 2.1867796
22 23 24 25 26 27 28
2.3402950 3.0690720 1.3446346 -0.1118490 0.1365966 -1.0192735 -1.6504919
29 30 31 32 33 34 35
-0.5379843 -2.3902736 -1.3638305 -3.3828945 -2.8800678 -2.7893295 -4.4230238
36 37 38 39 40 41 42
-3.5780075 -5.5078861 -6.3385707 -7.7364507 -4.3583950 -5.0866364 -4.5939771
43 44 45 46 47 48 49
-7.1132287 -5.5814236 -2.5700776 -2.7290657 -4.8720439 -0.3285659 1.0454413
50 51 52 53 54 55 56
4.7748348 8.9500509 5.7361274 6.2923542 6.6317892 6.2784023 5.2609552
57 58 59 60
2.1088253 -2.0575520 -0.9904444 -2.8290497
> postscript(file="/var/www/html/rcomp/tmp/68yju1258735091.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 0.9272595 NA
1 -2.8599596 0.9272595
2 -6.9322791 -2.8599596
3 -7.5553488 -6.9322791
4 -5.6828397 -7.5553488
5 -4.6198542 -5.6828397
6 -4.3625855 -4.6198542
7 -1.9298826 -4.3625855
8 1.1545405 -1.9298826
9 5.2356523 1.1545405
10 7.2164400 5.2356523
11 5.3909885 7.2164400
12 3.6470343 5.3909885
13 4.2870989 3.6470343
14 6.7379524 4.2870989
15 7.8281083 6.7379524
16 5.0151061 7.8281083
17 4.9723158 5.0151061
18 6.5612424 4.9723158
19 5.6332454 6.5612424
20 2.1867796 5.6332454
21 2.3402950 2.1867796
22 3.0690720 2.3402950
23 1.3446346 3.0690720
24 -0.1118490 1.3446346
25 0.1365966 -0.1118490
26 -1.0192735 0.1365966
27 -1.6504919 -1.0192735
28 -0.5379843 -1.6504919
29 -2.3902736 -0.5379843
30 -1.3638305 -2.3902736
31 -3.3828945 -1.3638305
32 -2.8800678 -3.3828945
33 -2.7893295 -2.8800678
34 -4.4230238 -2.7893295
35 -3.5780075 -4.4230238
36 -5.5078861 -3.5780075
37 -6.3385707 -5.5078861
38 -7.7364507 -6.3385707
39 -4.3583950 -7.7364507
40 -5.0866364 -4.3583950
41 -4.5939771 -5.0866364
42 -7.1132287 -4.5939771
43 -5.5814236 -7.1132287
44 -2.5700776 -5.5814236
45 -2.7290657 -2.5700776
46 -4.8720439 -2.7290657
47 -0.3285659 -4.8720439
48 1.0454413 -0.3285659
49 4.7748348 1.0454413
50 8.9500509 4.7748348
51 5.7361274 8.9500509
52 6.2923542 5.7361274
53 6.6317892 6.2923542
54 6.2784023 6.6317892
55 5.2609552 6.2784023
56 2.1088253 5.2609552
57 -2.0575520 2.1088253
58 -0.9904444 -2.0575520
59 -2.8290497 -0.9904444
60 NA -2.8290497
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.8599596 0.9272595
[2,] -6.9322791 -2.8599596
[3,] -7.5553488 -6.9322791
[4,] -5.6828397 -7.5553488
[5,] -4.6198542 -5.6828397
[6,] -4.3625855 -4.6198542
[7,] -1.9298826 -4.3625855
[8,] 1.1545405 -1.9298826
[9,] 5.2356523 1.1545405
[10,] 7.2164400 5.2356523
[11,] 5.3909885 7.2164400
[12,] 3.6470343 5.3909885
[13,] 4.2870989 3.6470343
[14,] 6.7379524 4.2870989
[15,] 7.8281083 6.7379524
[16,] 5.0151061 7.8281083
[17,] 4.9723158 5.0151061
[18,] 6.5612424 4.9723158
[19,] 5.6332454 6.5612424
[20,] 2.1867796 5.6332454
[21,] 2.3402950 2.1867796
[22,] 3.0690720 2.3402950
[23,] 1.3446346 3.0690720
[24,] -0.1118490 1.3446346
[25,] 0.1365966 -0.1118490
[26,] -1.0192735 0.1365966
[27,] -1.6504919 -1.0192735
[28,] -0.5379843 -1.6504919
[29,] -2.3902736 -0.5379843
[30,] -1.3638305 -2.3902736
[31,] -3.3828945 -1.3638305
[32,] -2.8800678 -3.3828945
[33,] -2.7893295 -2.8800678
[34,] -4.4230238 -2.7893295
[35,] -3.5780075 -4.4230238
[36,] -5.5078861 -3.5780075
[37,] -6.3385707 -5.5078861
[38,] -7.7364507 -6.3385707
[39,] -4.3583950 -7.7364507
[40,] -5.0866364 -4.3583950
[41,] -4.5939771 -5.0866364
[42,] -7.1132287 -4.5939771
[43,] -5.5814236 -7.1132287
[44,] -2.5700776 -5.5814236
[45,] -2.7290657 -2.5700776
[46,] -4.8720439 -2.7290657
[47,] -0.3285659 -4.8720439
[48,] 1.0454413 -0.3285659
[49,] 4.7748348 1.0454413
[50,] 8.9500509 4.7748348
[51,] 5.7361274 8.9500509
[52,] 6.2923542 5.7361274
[53,] 6.6317892 6.2923542
[54,] 6.2784023 6.6317892
[55,] 5.2609552 6.2784023
[56,] 2.1088253 5.2609552
[57,] -2.0575520 2.1088253
[58,] -0.9904444 -2.0575520
[59,] -2.8290497 -0.9904444
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.8599596 0.9272595
2 -6.9322791 -2.8599596
3 -7.5553488 -6.9322791
4 -5.6828397 -7.5553488
5 -4.6198542 -5.6828397
6 -4.3625855 -4.6198542
7 -1.9298826 -4.3625855
8 1.1545405 -1.9298826
9 5.2356523 1.1545405
10 7.2164400 5.2356523
11 5.3909885 7.2164400
12 3.6470343 5.3909885
13 4.2870989 3.6470343
14 6.7379524 4.2870989
15 7.8281083 6.7379524
16 5.0151061 7.8281083
17 4.9723158 5.0151061
18 6.5612424 4.9723158
19 5.6332454 6.5612424
20 2.1867796 5.6332454
21 2.3402950 2.1867796
22 3.0690720 2.3402950
23 1.3446346 3.0690720
24 -0.1118490 1.3446346
25 0.1365966 -0.1118490
26 -1.0192735 0.1365966
27 -1.6504919 -1.0192735
28 -0.5379843 -1.6504919
29 -2.3902736 -0.5379843
30 -1.3638305 -2.3902736
31 -3.3828945 -1.3638305
32 -2.8800678 -3.3828945
33 -2.7893295 -2.8800678
34 -4.4230238 -2.7893295
35 -3.5780075 -4.4230238
36 -5.5078861 -3.5780075
37 -6.3385707 -5.5078861
38 -7.7364507 -6.3385707
39 -4.3583950 -7.7364507
40 -5.0866364 -4.3583950
41 -4.5939771 -5.0866364
42 -7.1132287 -4.5939771
43 -5.5814236 -7.1132287
44 -2.5700776 -5.5814236
45 -2.7290657 -2.5700776
46 -4.8720439 -2.7290657
47 -0.3285659 -4.8720439
48 1.0454413 -0.3285659
49 4.7748348 1.0454413
50 8.9500509 4.7748348
51 5.7361274 8.9500509
52 6.2923542 5.7361274
53 6.6317892 6.2923542
54 6.2784023 6.6317892
55 5.2609552 6.2784023
56 2.1088253 5.2609552
57 -2.0575520 2.1088253
58 -0.9904444 -2.0575520
59 -2.8290497 -0.9904444
> 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/7nfj41258735091.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/848yh1258735091.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/9kvgo1258735091.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/10ahat1258735091.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/11afoh1258735091.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/12ug8c1258735091.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/1374op1258735091.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/14co3i1258735091.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/15cq6m1258735091.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/16yhmb1258735091.tab")
+ }
>
> system("convert tmp/14ygz1258735091.ps tmp/14ygz1258735091.png")
> system("convert tmp/2p79v1258735091.ps tmp/2p79v1258735091.png")
> system("convert tmp/3179d1258735091.ps tmp/3179d1258735091.png")
> system("convert tmp/4y3il1258735091.ps tmp/4y3il1258735091.png")
> system("convert tmp/5d8dx1258735091.ps tmp/5d8dx1258735091.png")
> system("convert tmp/68yju1258735091.ps tmp/68yju1258735091.png")
> system("convert tmp/7nfj41258735091.ps tmp/7nfj41258735091.png")
> system("convert tmp/848yh1258735091.ps tmp/848yh1258735091.png")
> system("convert tmp/9kvgo1258735091.ps tmp/9kvgo1258735091.png")
> system("convert tmp/10ahat1258735091.ps tmp/10ahat1258735091.png")
>
>
> proc.time()
user system elapsed
2.410 1.585 8.797