AST 431 L1
Monowar Hossain
2025-07-29
Econometrics 01
Lecture 01
library(haven)
#bd_birthweight <- read_sav("/home/mmhossain1/AST 431/bd_birthweight.sav")
bd_birthweight <- read_sav("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/bd_birthweight.sav")
#View(bd_birthweight)
library(haven)
bd_birthweight <- read_sav("Data/bd_birthweight.sav")
#View(bd_birthweight)
summary(bd_birthweight$bwt)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 886.2 2297.4 2589.3 2562.2 2826.3 4481.4y<-bd_birthweight$bwt
var(y)## [1] 152737.9sd(y)/mean(y)## [1] 0.1525305boxplot(y)
hist(y)
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsggplot(bd_birthweight) +
geom_boxplot(mapping = aes(x = bwt),fill = "blue" )
ggplot(bd_birthweight) +
geom_histogram(mapping = aes(x = bwt),fill ="green",color= "white")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
library(tidyverse)
ggplot(bd_birthweight) +
geom_density(aes(x = bwt),fill = "black")
## correlation (BIRTH WIGHT and BMI)
plot(bd_birthweight$bwt,bd_birthweight$v445)
cor(x = bd_birthweight$bwt,y = bd_birthweight$v445,use = "pairwise.complete.obs" )## [1] 0.03509002bd_birthweight%>%
select(bwt,v445)%>%
drop_na()%>%
cor()## bwt v445
## bwt 1.00000000 0.03509002
## v445 0.03509002 1.00000000 # correlation bet (bwt and number of child v218)
cor(bd_birthweight$bwt,bd_birthweight$v218)## [1] -0.02953598## boxplot
# v025= place od residence
ggplot(bd_birthweight) +
geom_boxplot(mapping = aes(y= bwt,group = v025))
low birth weight table
tab_f<-function(x,y){
r1<-table(x,y)
r2<-summary(r1)
r3<-prop.table(x = r1,margin = 2)*100
return(list(table=r1,summary = r2,Percentage =round(x = r3,digits = 2)))
}
#library(haven)
#bd_birthweight <- read_sav("C:/Users/training/Downloads/bd_birthweight.sav")
data<-bd_birthweight
# low = low birth weight
# v013 = age in 5-year groups
# v106 = highest educational level
# v701 = husband/partner's education level
# v190 = wealth index
# v704 = husband/partner's occupation
# v160 = toilet facilities shared with other households
# v025 = type of place of residence
# v024 = division
# v151 = sex of household head
# v504 = currently residing with husband/partner
# v209 = births in past year
tab_f(x = data$low,y = data$v013 )## $table
## y
## x 1 2 3 4 5 6 7
## 0 701 1628 1278 690 273 67 18
## 1 445 1086 886 533 214 46 21
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 11.098, df = 6, p-value = 0.0854
##
## $Percentage
## y
## x 1 2 3 4 5 6 7
## 0 61.17 59.99 59.06 56.42 56.06 59.29 46.15
## 1 38.83 40.01 40.94 43.58 43.94 40.71 53.85tab_f(x = data$low,y = data$v106 )## $table
## y
## x 0 1 2 3
## 0 656 1289 2189 521
## 1 577 917 1432 305
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 26.254, df = 3, p-value = 8.44e-06
##
## $Percentage
## y
## x 0 1 2 3
## 0 53.20 58.43 60.45 63.08
## 1 46.80 41.57 39.55 36.92tab_f(x = data$low,y = data$v701 )## $table
## y
## x 0 1 2 3
## 0 1091 1390 1448 725
## 1 917 987 912 414
##
## $summary
## Number of cases in table: 7884
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 33.96, df = 3, p-value = 2.018e-07
##
## $Percentage
## y
## x 0 1 2 3
## 0 54.33 58.48 61.36 63.65
## 1 45.67 41.52 38.64 36.35tab_f(x = data$low,y = data$v190 )## $table
## y
## x 1 2 3 4 5
## 0 971 892 885 969 938
## 1 766 611 631 633 590
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 12.279, df = 4, p-value = 0.01539
##
## $Percentage
## y
## x 1 2 3 4 5
## 0 55.90 59.35 58.38 60.49 61.39
## 1 44.10 40.65 41.62 39.51 38.61tab_f(x = data$low,y = data$v704 )## $table
## y
## x 11 12 13 14 15 16 21 22 23 31 41 51 52 61 62
## 0 17 589 401 106 4 7 523 8 374 1108 316 138 936 22 0
## 1 5 415 321 74 6 9 410 1 272 741 192 89 600 23 1
## y
## x 96 98
## 0 84 4
## 1 60 3
##
## $summary
## Number of cases in table: 7859
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 25.554, df = 16, p-value = 0.06064
## Chi-squared approximation may be incorrect
##
## $Percentage
## y
## x 11 12 13 14 15 16 21 22 23 31
## 0 77.27 58.67 55.54 58.89 40.00 43.75 56.06 88.89 57.89 59.92
## 1 22.73 41.33 44.46 41.11 60.00 56.25 43.94 11.11 42.11 40.08
## y
## x 41 51 52 61 62 96 98
## 0 62.20 60.79 60.94 48.89 0.00 58.33 57.14
## 1 37.80 39.21 39.06 51.11 100.00 41.67 42.86tab_f(x = data$low,y = data$v160 )## $table
## y
## x 0 1 7
## 0 2678 1466 404
## 1 1818 1059 267
##
## $summary
## Number of cases in table: 7692
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 1.871, df = 2, p-value = 0.3924
##
## $Percentage
## y
## x 0 1 7
## 0 59.56 58.06 60.21
## 1 40.44 41.94 39.79tab_f(x = data$low,y = data$v025 )## $table
## y
## x 1 2
## 0 1528 3127
## 1 960 2271
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 8.557, df = 1, p-value = 0.003442
##
## $Percentage
## y
## x 1 2
## 0 61.41 57.93
## 1 38.59 42.07tab_f(x = data$low,y = data$v024 )## $table
## y
## x 1 2 3 4 5 6 7
## 0 558 896 796 519 600 590 696
## 1 348 621 582 343 359 368 610
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 29.183, df = 6, p-value = 5.617e-05
##
## $Percentage
## y
## x 1 2 3 4 5 6 7
## 0 61.59 59.06 57.76 60.21 62.57 61.59 53.29
## 1 38.41 40.94 42.24 39.79 37.43 38.41 46.71tab_f(x = data$low,y = data$v151 )## $table
## y
## x 1 2
## 0 4229 426
## 1 2933 298
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 0.011759, df = 1, p-value = 0.9136
##
## $Percentage
## y
## x 1 2
## 0 59.05 58.84
## 1 40.95 41.16tab_f(x = data$low,y = data$v504 )## $table
## y
## x 1 2
## 0 3959 638
## 1 2759 420
##
## $summary
## Number of cases in table: 7776
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 0.7111, df = 1, p-value = 0.3991
##
## $Percentage
## y
## x 1 2
## 0 58.93 60.30
## 1 41.07 39.70tab_f(x = data$low,y = data$v209 )## $table
## y
## x 0 1 2
## 0 3338 1307 10
## 1 2416 803 12
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 11.547, df = 2, p-value = 0.003109
##
## $Percentage
## y
## x 0 1 2
## 0 58.01 61.94 45.45
## 1 41.99 38.06 54.55attr(data$low, "labels")## 2500-Hi 0-2500
## 0 1attr(data$v013, "labels")## 15-19 20-24 25-29 30-34 35-39 40-44 45-49
## 1 2 3 4 5 6 7attr(data$v106, "labels")## no education primary secondary higher
## 0 1 2 3attr(data$v701, "labels")## no education primary secondary higher don't know
## 0 1 2 3 8attr(data$v190, "labels")## poorest poorer middle richer richest
## 1 2 3 4 5attr(data$v704, "labels")## land owner
## 11
## farmer
## 12
## agricultural worker
## 13
## fisherman
## 14
## poultry raising, cattle raising
## 15
## home-based manufacturing (handicraft, food products)
## 16
## rickshaw driver, brick breaking, road building, construction worker, boatman, and earth work etc.
## 21
## domestic servant
## 22
## non-agricultural worker (factory worker, blue collar service)
## 23
## carpenter, mason, bus/taxi driver, construction supervisor, seamstresses/tailor, policeman, armed services, dai, commun
## 31
## doctor, lawyer, dentist, accountant, teacher, nurse, family welfare visitor, mid and high level services (government/pr
## 41
## big businessman
## 51
## small business/trader
## 52
## unemployed/student
## 61
## retired
## 62
## housewife
## 63
## beggar
## 64
## imam/religious leader
## 65
## others
## 96
## don't know
## 98attr(data$v160, "labels")## no yes not a dejure resident
## 0 1 7attr(data$v025, "labels")## urban rural
## 1 2attr(data$v024, "labels")## barisal chittagong dhaka khulna rajshahi rangpur sylhet
## 1 2 3 4 5 6 7attr(data$v151, "labels")## male female
## 1 2attr(data$v504, "labels")## living with her staying elsewhere
## 1 2attr(data$v209, "labels")## no births
## 0low birth we
model<- glm(formula = data$bwt ~ data$v013 + data$v106 + data$v701 + data$v190 + data$v704 +data$v160 + data$v025+data$v151+data$v504+data$v209,data = data)
summary(model)##
## Call:
## glm(formula = data$bwt ~ data$v013 + data$v106 + data$v701 +
## data$v190 + data$v704 + data$v160 + data$v025 + data$v151 +
## data$v504 + data$v209, data = data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2573.3944 32.3940 79.441 < 2e-16 ***
## data$v013 -3.5696 3.9057 -0.914 0.36078
## data$v106 9.1658 6.7521 1.357 0.17467
## data$v701 17.7185 5.9413 2.982 0.00287 **
## data$v190 1.0497 4.1816 0.251 0.80180
## data$v704 0.1101 0.2812 0.391 0.69552
## data$v160 -0.3352 2.3449 -0.143 0.88634
## data$v025 -15.1288 10.7293 -1.410 0.15857
## data$v151 -14.9621 18.8917 -0.792 0.42839
## data$v504 -6.9974 15.7139 -0.445 0.65612
## data$v209 14.7431 10.0626 1.465 0.14293
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 151317.8)
##
## Null deviance: 1148249085 on 7555 degrees of freedom
## Residual deviance: 1141692425 on 7545 degrees of freedom
## (330 observations deleted due to missingness)
## AIC: 111577
##
## Number of Fisher Scoring iterations: 2# Tidy output: Coefficients (β), SE, p-values, and 95% CI for OR
tidy_model <- broom::tidy(model, conf.int = T)
library(tidyverse)
# Display the model output in a readable format
tidy_model %>%
select(term, estimate, std.error, statistic, p.value, conf.low, conf.high) %>%
mutate(OR = estimate) ## # A tibble: 11 × 8
## term estimate std.error statistic p.value conf.low conf.high OR
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 2573. 32.4 79.4 0 2510. 2637. 2573.
## 2 data$v013 -3.57 3.91 -0.914 0.361 -11.2 4.09 -3.57
## 3 data$v106 9.17 6.75 1.36 0.175 -4.07 22.4 9.17
## 4 data$v701 17.7 5.94 2.98 0.00287 6.07 29.4 17.7
## 5 data$v190 1.05 4.18 0.251 0.802 -7.15 9.25 1.05
## 6 data$v704 0.110 0.281 0.391 0.696 -0.441 0.661 0.110
## 7 data$v160 -0.335 2.34 -0.143 0.886 -4.93 4.26 -0.335
## 8 data$v025 -15.1 10.7 -1.41 0.159 -36.2 5.90 -15.1
## 9 data$v151 -15.0 18.9 -0.792 0.428 -52.0 22.1 -15.0
## 10 data$v504 -7.00 15.7 -0.445 0.656 -37.8 23.8 -7.00
## 11 data$v209 14.7 10.1 1.47 0.143 -4.98 34.5 14.7Lecture 02
class work
cross table & association(chi_squrae)
tab_f<-function(x,y){
r1<-table(x,y)
r2<-summary(r1)
r3<-prop.table(x = r1,margin = 2)*100
return(list(table=r1,summary = r2,Percentage =round(x = r3,digits = 2)))
}
#library(haven)
#bd_birthweight <- read_sav("C:/Users/training/Downloads/bd_birthweight.sav")
data<-bd_birthweight
# low = low birth weight
# v013 = age in 5-year groups
# v106 = highest educational level
# v701 = husband/partner's education level
# v190 = wealth index
# v704 = husband/partner's occupation
# v160 = toilet facilities shared with other households
# v025 = type of place of residence
# v024 = division
# v151 = sex of household head
# v504 = currently residing with husband/partner
# v209 = births in past year
tab_f(x = data$low,y = data$v013 )## $table
## y
## x 1 2 3 4 5 6 7
## 0 701 1628 1278 690 273 67 18
## 1 445 1086 886 533 214 46 21
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 11.098, df = 6, p-value = 0.0854
##
## $Percentage
## y
## x 1 2 3 4 5 6 7
## 0 61.17 59.99 59.06 56.42 56.06 59.29 46.15
## 1 38.83 40.01 40.94 43.58 43.94 40.71 53.85tab_f(x = data$low,y = data$v106 )## $table
## y
## x 0 1 2 3
## 0 656 1289 2189 521
## 1 577 917 1432 305
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 26.254, df = 3, p-value = 8.44e-06
##
## $Percentage
## y
## x 0 1 2 3
## 0 53.20 58.43 60.45 63.08
## 1 46.80 41.57 39.55 36.92tab_f(x = data$low,y = data$v701 )## $table
## y
## x 0 1 2 3
## 0 1091 1390 1448 725
## 1 917 987 912 414
##
## $summary
## Number of cases in table: 7884
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 33.96, df = 3, p-value = 2.018e-07
##
## $Percentage
## y
## x 0 1 2 3
## 0 54.33 58.48 61.36 63.65
## 1 45.67 41.52 38.64 36.35tab_f(x = data$low,y = data$v190 )## $table
## y
## x 1 2 3 4 5
## 0 971 892 885 969 938
## 1 766 611 631 633 590
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 12.279, df = 4, p-value = 0.01539
##
## $Percentage
## y
## x 1 2 3 4 5
## 0 55.90 59.35 58.38 60.49 61.39
## 1 44.10 40.65 41.62 39.51 38.61tab_f(x = data$low,y = data$v704 )## $table
## y
## x 11 12 13 14 15 16 21 22 23 31 41 51 52 61 62
## 0 17 589 401 106 4 7 523 8 374 1108 316 138 936 22 0
## 1 5 415 321 74 6 9 410 1 272 741 192 89 600 23 1
## y
## x 96 98
## 0 84 4
## 1 60 3
##
## $summary
## Number of cases in table: 7859
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 25.554, df = 16, p-value = 0.06064
## Chi-squared approximation may be incorrect
##
## $Percentage
## y
## x 11 12 13 14 15 16 21 22 23 31
## 0 77.27 58.67 55.54 58.89 40.00 43.75 56.06 88.89 57.89 59.92
## 1 22.73 41.33 44.46 41.11 60.00 56.25 43.94 11.11 42.11 40.08
## y
## x 41 51 52 61 62 96 98
## 0 62.20 60.79 60.94 48.89 0.00 58.33 57.14
## 1 37.80 39.21 39.06 51.11 100.00 41.67 42.86tab_f(x = data$low,y = data$v160 )## $table
## y
## x 0 1 7
## 0 2678 1466 404
## 1 1818 1059 267
##
## $summary
## Number of cases in table: 7692
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 1.871, df = 2, p-value = 0.3924
##
## $Percentage
## y
## x 0 1 7
## 0 59.56 58.06 60.21
## 1 40.44 41.94 39.79tab_f(x = data$low,y = data$v025 )## $table
## y
## x 1 2
## 0 1528 3127
## 1 960 2271
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 8.557, df = 1, p-value = 0.003442
##
## $Percentage
## y
## x 1 2
## 0 61.41 57.93
## 1 38.59 42.07tab_f(x = data$low,y = data$v024 )## $table
## y
## x 1 2 3 4 5 6 7
## 0 558 896 796 519 600 590 696
## 1 348 621 582 343 359 368 610
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 29.183, df = 6, p-value = 5.617e-05
##
## $Percentage
## y
## x 1 2 3 4 5 6 7
## 0 61.59 59.06 57.76 60.21 62.57 61.59 53.29
## 1 38.41 40.94 42.24 39.79 37.43 38.41 46.71tab_f(x = data$low,y = data$v151 )## $table
## y
## x 1 2
## 0 4229 426
## 1 2933 298
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 0.011759, df = 1, p-value = 0.9136
##
## $Percentage
## y
## x 1 2
## 0 59.05 58.84
## 1 40.95 41.16tab_f(x = data$low,y = data$v504 )## $table
## y
## x 1 2
## 0 3959 638
## 1 2759 420
##
## $summary
## Number of cases in table: 7776
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 0.7111, df = 1, p-value = 0.3991
##
## $Percentage
## y
## x 1 2
## 0 58.93 60.30
## 1 41.07 39.70tab_f(x = data$low,y = data$v209 )## $table
## y
## x 0 1 2
## 0 3338 1307 10
## 1 2416 803 12
##
## $summary
## Number of cases in table: 7886
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 11.547, df = 2, p-value = 0.003109
##
## $Percentage
## y
## x 0 1 2
## 0 58.01 61.94 45.45
## 1 41.99 38.06 54.55attr(data$low, "labels")## 2500-Hi 0-2500
## 0 1attr(data$v013, "labels")## 15-19 20-24 25-29 30-34 35-39 40-44 45-49
## 1 2 3 4 5 6 7attr(data$v106, "labels")## no education primary secondary higher
## 0 1 2 3attr(data$v701, "labels")## no education primary secondary higher don't know
## 0 1 2 3 8attr(data$v190, "labels")## poorest poorer middle richer richest
## 1 2 3 4 5attr(data$v704, "labels")## land owner
## 11
## farmer
## 12
## agricultural worker
## 13
## fisherman
## 14
## poultry raising, cattle raising
## 15
## home-based manufacturing (handicraft, food products)
## 16
## rickshaw driver, brick breaking, road building, construction worker, boatman, and earth work etc.
## 21
## domestic servant
## 22
## non-agricultural worker (factory worker, blue collar service)
## 23
## carpenter, mason, bus/taxi driver, construction supervisor, seamstresses/tailor, policeman, armed services, dai, commun
## 31
## doctor, lawyer, dentist, accountant, teacher, nurse, family welfare visitor, mid and high level services (government/pr
## 41
## big businessman
## 51
## small business/trader
## 52
## unemployed/student
## 61
## retired
## 62
## housewife
## 63
## beggar
## 64
## imam/religious leader
## 65
## others
## 96
## don't know
## 98attr(data$v160, "labels")## no yes not a dejure resident
## 0 1 7attr(data$v025, "labels")## urban rural
## 1 2attr(data$v024, "labels")## barisal chittagong dhaka khulna rajshahi rangpur sylhet
## 1 2 3 4 5 6 7attr(data$v151, "labels")## male female
## 1 2attr(data$v504, "labels")## living with her staying elsewhere
## 1 2attr(data$v209, "labels")## no births
## 0linear regression (cont .. var = bwt)
# v012 = age in 5-year groups
# v025 = type of place of residence
# v106 = highest educational level
# v190 = wealth index
# v701 = husband/partner's education level
model<-lm(formula = bwt~ v012 + as.factor(v025)+ as.factor(v106) + as.factor(v190) + as.factor(v701),data = data)
model##
## Call:
## lm(formula = bwt ~ v012 + as.factor(v025) + as.factor(v106) +
## as.factor(v190) + as.factor(v701), data = data)
##
## Coefficients:
## (Intercept) v012 as.factor(v025)2 as.factor(v106)1
## 2537.101 -1.032 -14.510 35.496
## as.factor(v106)2 as.factor(v106)3 as.factor(v190)2 as.factor(v190)3
## 31.448 37.382 6.900 -11.273
## as.factor(v190)4 as.factor(v190)5 as.factor(v701)1 as.factor(v701)2
## -3.107 12.177 35.541 43.070
## as.factor(v701)3
## 59.966summary(model) ##
## Call:
## lm(formula = bwt ~ v012 + as.factor(v025) + as.factor(v106) +
## as.factor(v190) + as.factor(v701), data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1676.92 -260.45 26.09 261.12 1921.66
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2537.1007 27.9574 90.749 < 2e-16 ***
## v012 -1.0325 0.7971 -1.295 0.19528
## as.factor(v025)2 -14.5104 10.7181 -1.354 0.17584
## as.factor(v106)1 35.4962 14.5825 2.434 0.01495 *
## as.factor(v106)2 31.4475 15.4153 2.040 0.04138 *
## as.factor(v106)3 37.3819 22.3945 1.669 0.09511 .
## as.factor(v190)2 6.9001 13.9534 0.495 0.62096
## as.factor(v190)3 -11.2725 14.4592 -0.780 0.43564
## as.factor(v190)4 -3.1070 15.1209 -0.205 0.83720
## as.factor(v190)5 12.1765 17.4774 0.697 0.48601
## as.factor(v701)1 35.5407 12.5005 2.843 0.00448 **
## as.factor(v701)2 43.0696 13.9908 3.078 0.00209 **
## as.factor(v701)3 59.9656 19.2051 3.122 0.00180 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 389.8 on 7871 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.007179, Adjusted R-squared: 0.005665
## F-statistic: 4.743 on 12 and 7871 DF, p-value: 8.791e-08Time Series 01
Lecture 03
#births <- read.table("~/AST 431/births.txt", quote="\"", comment.char="")
births <- read.table("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/births.txt", quote="\"", comment.char="")
#View(births)
#Kings <- read.table("~/AST 431/Kings.txt", quote="\"", comment.char="")
Kings <- read.table("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/Kings.txt", quote="\"", comment.char="")
#View(Kings)
kings_time_series<-ts(Kings)
plot(Kings)
plot(kings_time_series)
births_time_series<-ts(births,frequency = 12,start = c(1946,1))
births_time_series## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
## 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
## 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
## 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
## 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
## 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
## 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
## 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
## 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
## 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
## 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
## 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
## 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
## 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
## Nov Dec
## 1946 21.672 21.870
## 1947 21.759 22.073
## 1948 21.059 21.573
## 1949 21.519 22.025
## 1950 22.084 22.991
## 1951 22.964 23.981
## 1952 23.162 24.707
## 1953 25.246 25.180
## 1954 24.712 25.688
## 1955 25.693 26.881
## 1956 26.291 26.987
## 1957 26.634 27.735
## 1958 25.912 26.619
## 1959 26.992 27.897plot(births)
plot(births_time_series)
decompose(births_time_series)## $x
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
## 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
## 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
## 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
## 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
## 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
## 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
## 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
## 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
## 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
## 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
## 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
## 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
## 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
## Nov Dec
## 1946 21.672 21.870
## 1947 21.759 22.073
## 1948 21.059 21.573
## 1949 21.519 22.025
## 1950 22.084 22.991
## 1951 22.964 23.981
## 1952 23.162 24.707
## 1953 25.246 25.180
## 1954 24.712 25.688
## 1955 25.693 26.881
## 1956 26.291 26.987
## 1957 26.634 27.735
## 1958 25.912 26.619
## 1959 26.992 27.897
##
## $seasonal
## Jan Feb Mar Apr May Jun
## 1946 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1947 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1948 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1949 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1950 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1951 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1952 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1953 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1954 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1955 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1956 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1957 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1958 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## 1959 -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## Jul Aug Sep Oct Nov Dec
## 1946 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1947 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1948 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1949 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1950 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1951 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1952 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1953 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1954 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1955 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1956 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1957 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1958 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
## 1959 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
##
## $trend
## Jan Feb Mar Apr May Jun Jul Aug
## 1946 NA NA NA NA NA NA 23.98433 23.66213
## 1947 22.35350 22.30871 22.30258 22.29479 22.29354 22.30562 22.33483 22.31167
## 1948 22.43038 22.43667 22.38721 22.35242 22.32458 22.27458 22.23754 22.21988
## 1949 22.06375 22.08033 22.13317 22.16604 22.17542 22.21342 22.27625 22.35750
## 1950 23.21663 23.26967 23.33492 23.42679 23.50638 23.57017 23.63888 23.75713
## 1951 24.00083 24.12350 24.20917 24.28208 24.35450 24.43242 24.49496 24.48379
## 1952 24.27204 24.27300 24.28942 24.30129 24.31325 24.35175 24.40558 24.44475
## 1953 24.78646 24.84992 24.92692 25.02362 25.16308 25.26963 25.30154 25.34125
## 1954 25.92446 25.92317 25.92967 25.92137 25.89567 25.89458 25.92963 25.98246
## 1955 25.64612 25.78679 25.93192 26.06388 26.16329 26.25388 26.35471 26.40496
## 1956 27.21104 27.21900 27.20700 27.26925 27.35050 27.37983 27.39975 27.44150
## 1957 27.44221 27.40283 27.44300 27.45717 27.44429 27.48975 27.54354 27.56933
## 1958 27.68642 27.76067 27.75963 27.71037 27.65783 27.58125 27.49075 27.46183
## 1959 26.96858 27.00512 27.09250 27.17263 27.26208 27.36033 NA NA
## Sep Oct Nov Dec
## 1946 23.42333 23.16112 22.86425 22.54521
## 1947 22.26279 22.25796 22.27767 22.35400
## 1948 22.16983 22.07721 22.01396 22.02604
## 1949 22.48862 22.70992 22.98563 23.16346
## 1950 23.86354 23.89533 23.87342 23.88150
## 1951 24.43879 24.36829 24.29192 24.27642
## 1952 24.49325 24.58517 24.70429 24.76017
## 1953 25.42779 25.57588 25.73904 25.87513
## 1954 26.01054 25.88617 25.67087 25.57312
## 1955 26.45379 26.64933 26.95183 27.14683
## 1956 27.45229 27.43354 27.44488 27.46996
## 1957 27.63167 27.67804 27.62579 27.61212
## 1958 27.42262 27.34175 27.25129 27.08558
## 1959 NA NA NA NA
##
## $random
## Jan Feb Mar Apr May
## 1946 NA NA NA NA NA
## 1947 -0.237305288 0.863252404 0.543893429 0.175887019 -0.793193109
## 1948 0.183819712 -0.318705929 0.340268429 0.121262019 -0.354234776
## 1949 0.161444712 0.002627404 -0.571689904 -0.749362981 -0.666068109
## 1950 0.064569712 -0.292705929 0.479560096 1.047887019 1.561973558
## 1951 -0.036638622 1.008460737 0.004310096 0.556595353 -0.176151442
## 1952 0.203153045 0.079960737 -0.376939904 -0.853612981 -0.576901442
## 1953 0.254736378 -0.122955929 -0.224439904 -0.159946314 0.016265224
## 1954 -0.590263622 -0.536205929 0.189810096 1.079303686 1.062681891
## 1955 0.021069712 0.535169071 -0.073439904 -1.787196314 -1.647943109
## 1956 -0.316846955 -0.918039263 -0.155523237 0.507428686 0.924848558
## 1957 -0.176013622 -0.471872596 -0.762523237 0.240512019 1.182056891
## 1958 0.122778045 -0.753705929 0.340851763 -0.319696314 0.021515224
## 1959 -0.215388622 0.363835737 -0.295023237 -0.419946314 -1.115734776
## Jun Jul Aug Sep Oct
## 1946 NA -0.963379006 -0.925718750 -0.939949519 -0.709369391
## 1947 -1.391369391 -0.311879006 0.347739583 0.150592147 0.076797276
## 1948 0.001672276 0.256412660 0.119531250 -0.623449519 0.289547276
## 1949 0.813838942 0.371704327 0.225906250 0.081758814 -0.578161058
## 1950 0.166088942 -0.423920673 -0.467718750 -0.433157853 -0.418577724
## 1951 0.387838942 0.499995994 -0.030385417 -0.116407853 -0.033536058
## 1952 0.538505609 0.414370994 0.206656250 0.025133814 -0.161411058
## 1953 -0.481369391 0.251412660 0.100156250 0.148592147 0.110880609
## 1954 0.380672276 -0.679670673 -0.269052083 -0.550157853 -0.282411058
## 1955 0.118380609 0.550245994 1.029447917 0.768592147 0.359422276
## 1956 -0.087577724 0.126204327 -0.437093750 -0.087907853 0.927213942
## 1957 0.053505609 -0.934587340 -0.592927083 0.724717147 0.030713942
## 1958 0.581005609 0.282204327 0.132572917 0.290758814 -0.171994391
## 1959 -1.642077724 NA NA NA NA
## Nov Dec
## 1946 -0.082484776 -0.298388622
## 1947 0.591098558 0.095819712
## 1948 0.154806891 -0.076221955
## 1949 -0.356859776 -0.761638622
## 1950 -0.679651442 -0.513680288
## 1951 -0.218151442 0.081403045
## 1952 -0.432526442 0.323653045
## 1953 0.616723558 -0.318305288
## 1954 0.150890224 0.491694712
## 1955 -0.149068109 0.110986378
## 1956 -0.044109776 -0.106138622
## 1957 0.117973558 0.499694712
## 1958 -0.229526442 -0.089763622
## 1959 NA NA
##
## $figure
## [1] -0.6771947 -2.0829607 0.8625232 -0.8016787 0.2516514 -0.1532556
## [7] 1.4560457 1.1645938 0.6916162 0.7752444 -1.1097652 -0.3768197
##
## $type
## [1] "additive"
##
## attr(,"class")
## [1] "decomposed.ts"plot(decompose(births_time_series))
Lecture 04
library(readr)
#Rainfall_dhaka <- read_csv("/home/mmhossain1/AST 431/Rainfall_dhaka.csv")
Rainfall_dhaka <- read_csv("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/Rainfall_dhaka2.csv")## Rows: 526 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Stati, 14, 15, 31
## dbl (30): Year, Month, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 17, 18...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#View(Rainfall_dhaka)
library(ggplot2)
head(Rainfall_dhaka)## # A tibble: 6 × 34
## Stati Year Month `1` `2` `3` `4` `5` `6` `7` `8` `9` `10`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Dhaka 1976 1 0 0 0 0 0 0 0 0 0 0
## 2 Dhaka 1976 2 0 1 0 0 2 3 1 0 0 0
## 3 Dhaka 1976 3 0 0 0 7 0 0 0 34 0 0
## 4 Dhaka 1976 4 0 0 0 1 0 0 0 0 0 0
## 5 Dhaka 1976 5 0 0 0 0 1 12 1 24 24 76
## 6 Dhaka 1976 6 0 43 0 14 0 1 0 163 100 66
## # ℹ 21 more variables: `11` <dbl>, `12` <dbl>, `13` <dbl>, `14` <chr>,
## # `15` <chr>, `16` <dbl>, `17` <dbl>, `18` <dbl>, `19` <dbl>, `20` <dbl>,
## # `21` <dbl>, `22` <dbl>, `23` <dbl>, `24` <dbl>, `25` <dbl>, `26` <dbl>,
## # `27` <dbl>, `28` <dbl>, `29` <dbl>, `30` <dbl>, `31` <chr>#basic line plot
ggplot(data = economics,mapping = aes(x = date ,y=pop ))+
geom_line( color ="#00AFBB")
#ggplot(data = Rainfall_dhaka,mapping = aes(x = Year ,y= mean))
library(tidyverse, lib.loc = "/usr/lib/R/site-library")
library(dplyr)data cleaning
#View(Rainfall_dhaka)
# Replace values in the column `14` where the condition holds (e.g., values equal to some number)
Rainfall_dhaka$`14`[Rainfall_dhaka$`14` == "***"] <- NA
Rainfall_dhaka$`15`[Rainfall_dhaka$`15` == "***"] <- NA
#Rainfall_dhaka$`31`[Rainfall_dhaka$`31` == "***"] <- NA
Rainfall_dhaka$`14`<-as.double(Rainfall_dhaka$`14`)
Rainfall_dhaka$`15`<-as.double(Rainfall_dhaka$`15`)
Rainfall_dhaka$`31`<-as.double(Rainfall_dhaka$`31`)## Warning: NAs introduced by coercionRainfall_dhaka$mean<-rowMeans(Rainfall_dhaka[,4:34],na.rm = T)
Rainfall_dhaka$mean## [1] 0.00000000 0.24137931 3.77419355 1.13333333 14.80645161 20.90000000
## [7] 11.16129032 11.64516129 5.50000000 3.67741935 0.26666667 0.00000000
## [13] 0.00000000 2.35714286 2.29032258 8.50000000 12.29032258 8.40000000
## [19] 9.87096774 2.96774194 4.36666667 8.80645161 0.33333333 0.77419355
## [25] 0.00000000 0.71428571 0.58064516 6.46666667 14.64516129 17.63333333
## [31] 10.32258065 13.74193548 6.40000000 3.16129032 0.00000000 0.00000000
## [37] 0.09677419 0.46428571 0.19354839 0.56666667 3.67741935 8.60000000
## [43] 8.61290323 16.93548387 12.73333333 4.70967742 1.83333333 1.64516129
## [49] 0.09677419 1.10344828 1.74193548 4.90000000 13.35483871 10.76666667
## [55] 12.25806452 8.67741935 9.86666667 9.67741935 0.00000000 0.00000000
## [61] 0.32258065 1.50000000 3.51612903 9.13333333 8.77419355 5.60000000
## [67] 11.48387097 6.06451613 10.66666667 2.64516129 0.30000000 1.12903226
## [73] 0.00000000 0.53571429 2.61290323 3.46666667 4.96774194 17.13333333
## [79] 4.38709677 11.16129032 8.60000000 4.70967742 1.70000000 0.00000000
## [85] 0.46666667 2.17857143 4.45161290 10.60000000 11.22580645 10.00000000
## [91] 5.77419355 14.09677419 10.73333333 8.16129032 0.00000000 0.58064516
## [97] 0.41935484 0.03448276 0.16129032 4.13333333 22.80645161 21.23333333
## [103] 22.38709677 10.03225806 15.93333333 1.87096774 0.00000000 0.00000000
## [109] 0.25806452 0.03571429 6.29032258 5.86666667 9.67741935 13.30000000
## [115] 8.45161290 10.22580645 10.20000000 2.54838710 0.00000000 0.32258065
## [121] 0.70967742 0.00000000 0.74193548 8.23333333 6.16129032 10.13333333
## [127] 14.29032258 5.51612903 22.90000000 7.64516129 5.73333333 0.09677419
## [133] 0.12903226 0.00000000 1.06451613 7.66666667 3.51612903 10.53333333
## [139] 16.96774194 14.90322581 12.10000000 3.35483871 0.23333333 1.06451613
## [145] 0.00000000 1.51724138 2.38709677 9.40000000 16.54838710 19.33333333
## [151] 8.22580645 5.45161290 6.53333333 6.87096774 5.10000000 0.09677419
## [157] 0.00000000 1.14285714 0.00000000 2.83333333 7.35483871 10.63333333
## [163] 11.19354839 1.90322581 10.16666667 7.74193548 0.00000000 0.38709677
## [169] 0.00000000 1.28571429 4.87096774 5.13333333 6.51612903 7.63333333
## [175] 18.29032258 7.32258065 8.23333333 5.83870968 3.43333333 0.19354839
## [181] 0.87096774 0.28571429 1.48387097 1.76666667 17.06451613 10.66666667
## [187] 10.25806452 11.12903226 23.06666667 12.64516129 0.46666667 3.41935484
## [193] 0.03225806 1.62068966 0.00000000 0.83333333 4.93548387 4.40000000
## [199] 12.45161290 5.87096774 5.26666667 2.67741935 0.06666667 0.00000000
## [205] 0.00000000 1.85714286 2.83870968 3.76666667 17.93548387 16.80000000
## [211] 13.58064516 13.93548387 13.90000000 7.00000000 0.63333333 0.00000000
## [217] 0.41935484 1.92857143 3.70967742 6.70000000 8.19354839 8.86666667
## [223] 4.93548387 7.93548387 5.63333333 1.77419355 0.46666667 0.00000000
## [229] 0.25806452 1.10714286 0.00000000 2.93333333 8.51612903 7.90000000
## [235] 11.41935484 11.61290323 6.83333333 2.93548387 3.73333333 0.03225806
## [241] 0.00000000 0.72413793 1.74193548 6.63333333 6.70967742 11.43333333
## [247] 8.29032258 11.64516129 8.13333333 11.51612903 0.00000000 0.00000000
## [253] 0.06451613 0.25000000 2.64516129 4.43333333 4.87096774 8.30000000
## [259] 17.70967742 7.41935484 14.66666667 0.96774194 0.03333333 0.70967742
## [265] 1.58064516 0.14285714 2.67741935 5.93333333 13.06451613 2.96666667
## [271] 16.80645161 17.80645161 8.20000000 3.22580645 2.76666667 0.00000000
## [277] 0.00000000 0.00000000 0.00000000 0.70000000 13.80645161 11.60000000
## [283] 17.83870968 9.09677419 12.03333333 11.87096774 0.43333333 0.00000000
## [289] 0.41935484 1.51724138 5.54838710 6.30000000 19.61290323 5.50000000
## [295] 6.35483871 11.58064516 7.20000000 8.96774194 0.00000000 0.00000000
## [301] 0.00000000 0.03571429 1.06451613 1.53333333 12.96774194 12.86666667
## [307] 6.51612903 6.61290323 6.96666667 5.70967742 0.60000000 0.00000000
## [313] 0.70967742 0.14285714 1.64516129 3.70000000 8.77419355 12.43333333
## [319] 14.38709677 8.77419355 5.20000000 1.67741935 1.20000000 0.00000000
## [325] 0.00000000 0.89285714 3.09677419 4.10000000 4.51612903 15.76666667
## [331] 6.16129032 6.51612903 8.80000000 4.32258065 0.00000000 1.45161290
## [337] 0.00000000 0.00000000 0.29032258 5.56666667 5.22580645 15.86666667
## [343] 9.51612903 6.16129032 27.96666667 6.70967742 0.00000000 0.00000000
## [349] 0.03225806 0.10714286 5.00000000 3.03333333 9.38709677 8.63333333
## [355] 17.48387097 11.64516129 17.13333333 13.45161290 0.10000000 0.00000000
## [361] 0.00000000 0.00000000 0.00000000 6.03333333 5.96774194 10.86666667
## [367] 10.67741935 5.38709677 22.10000000 1.96774194 0.16666667 0.00000000
## [373] 0.00000000 1.07142857 0.35483871 5.43333333 5.96774194 20.93333333
## [379] 24.29032258 16.29032258 5.96666667 10.32258065 3.70000000 0.00000000
## [385] 0.74193548 1.93103448 1.45161290 3.03333333 6.61290323 13.63333333
## [391] 18.16129032 10.29032258 9.30000000 7.32258065 0.00000000 0.00000000
## [397] 0.03225806 0.03571429 1.38709677 0.46666667 5.41935484 5.66666667
## [403] 21.80645161 15.54838710 9.93333333 2.38709677 0.13333333 0.00000000
## [409] 0.00000000 1.71428571 0.70967742 1.23333333 5.70967742 10.26666667
## [415] 5.38709677 10.96774194 5.63333333 5.61290323 0.00000000 2.61290323
## [421] 0.00000000 0.00000000 0.64516129 4.10000000 7.58064516 10.46666667
## [427] 11.48387097 13.19354839 6.90000000 3.61290323 0.00000000 0.00000000
## [433] 0.32258065 0.03448276 1.19354839 8.96666667 4.41935484 5.83333333
## [439] 7.29032258 9.09677419 2.70000000 1.22580645 2.26666667 0.16129032
## [445] 0.00000000 0.28571429 0.83870968 1.06666667 12.19354839 10.83333333
## [451] 9.74193548 6.83870968 4.60000000 4.22580645 0.00000000 0.12903226
## [457] 0.00000000 0.42857143 0.32258065 2.66666667 4.74193548 11.40000000
## [463] 6.83870968 12.61290323 5.20000000 1.58064516 0.00000000 0.00000000
## [469] 0.09677419 0.60714286 0.12903226 5.53333333 5.96774194 12.50000000
## [475] 20.09677419 12.74193548 11.53333333 1.64516129 0.00000000 0.03225806
## [481] 0.09677419 0.44827586 1.77419355 1.83333333 6.83870968 7.06666667
## [487] 13.06451613 5.51612903 4.60000000 2.45161290 0.83333333 0.00000000
## [493] 0.00000000 0.07142857 3.22580645 7.60000000 6.06451613 13.80000000
## [499] 18.83870968 17.54838710 12.70000000 13.29032258 0.20000000 1.06451613
## [505] 0.00000000 0.71428571 0.09677419 10.30000000 12.64516129 12.20000000
## [511] 11.41935484 4.54838710 2.53333333 1.45161290 0.43333333 0.41935484
## [517] 0.03225806 4.10714286 1.25806452 7.06666667 7.45161290 8.06666667
## [523] 12.35483871 7.19354839 5.36666667 6.06451613#View(Rainfall_dhaka)
# Replace values in the column `14` where the condition holds (e.g., values equal to some number)
Rainfall_dhaka$`14`[Rainfall_dhaka$`14` == "***"] <- NA
Rainfall_dhaka$`15`[Rainfall_dhaka$`15` == "***"] <- NA
#Rainfall_dhaka$`31`[Rainfall_dhaka$`31` == "***"] <- NA
Rainfall_dhaka$`14`<-as.double(Rainfall_dhaka$`14`)
Rainfall_dhaka$`15`<-as.double(Rainfall_dhaka$`15`)
Rainfall_dhaka$`31`<-as.double(Rainfall_dhaka$`31`)
Rainfall_dhaka %>%
select(4:34) %>%
mutate(mean2= rowMeans(Rainfall_dhaka[,4:34],na.rm=T))## # A tibble: 526 × 32
## `1` `2` `3` `4` `5` `6` `7` `8` `9` `10` `11` `12` `13`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 1 0 0 2 3 1 0 0 0 0 0 0
## 3 0 0 0 7 0 0 0 34 0 0 0 0 0
## 4 0 0 0 1 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 1 12 1 24 24 76 6 17 5
## 6 0 43 0 14 0 1 0 163 100 66 87 20 0
## 7 0 33 10 11 4 0 0 103 0 52 28 28 0
## 8 0 6 12 20 11 9 4 6 0 0 0 0 6
## 9 13 0 0 0 0 0 0 0 0 3 0 3 0
## 10 17 18 1 6 5 0 27 0 0 0 0 0 0
## # ℹ 516 more rows
## # ℹ 19 more variables: `14` <dbl>, `15` <dbl>, `16` <dbl>, `17` <dbl>,
## # `18` <dbl>, `19` <dbl>, `20` <dbl>, `21` <dbl>, `22` <dbl>, `23` <dbl>,
## # `24` <dbl>, `25` <dbl>, `26` <dbl>, `27` <dbl>, `28` <dbl>, `29` <dbl>,
## # `30` <dbl>, `31` <dbl>, mean2 <dbl>Lecture 06
library(readr)
#rainfall <- read_csv("/home/mmhossain1/AST 431/Rainfall_dhaka.csv")
rainfall <- read_csv("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/Rainfall_dhaka2.csv")## Rows: 526 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Stati, 14, 15, 31
## dbl (30): Year, Month, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 17, 18...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.rainfall <- rainfall %>% select(!(Stati))
str(rainfall)## tibble [526 × 33] (S3: tbl_df/tbl/data.frame)
## $ Year : num [1:526] 1976 1976 1976 1976 1976 ...
## $ Month: num [1:526] 1 2 3 4 5 6 7 8 9 10 ...
## $ 1 : num [1:526] 0 0 0 0 0 0 0 0 13 17 ...
## $ 2 : num [1:526] 0 1 0 0 0 43 33 6 0 18 ...
## $ 3 : num [1:526] 0 0 0 0 0 0 10 12 0 1 ...
## $ 4 : num [1:526] 0 0 7 1 0 14 11 20 0 6 ...
## $ 5 : num [1:526] 0 2 0 0 1 0 4 11 0 5 ...
## $ 6 : num [1:526] 0 3 0 0 12 1 0 9 0 0 ...
## $ 7 : num [1:526] 0 1 0 0 1 0 0 4 0 27 ...
## $ 8 : num [1:526] 0 0 34 0 24 163 103 6 0 0 ...
## $ 9 : num [1:526] 0 0 0 0 24 100 0 0 0 0 ...
## $ 10 : num [1:526] 0 0 0 0 76 66 52 0 3 0 ...
## $ 11 : num [1:526] 0 0 0 0 6 87 28 0 0 0 ...
## $ 12 : num [1:526] 0 0 0 0 17 20 28 0 3 0 ...
## $ 13 : num [1:526] 0 0 0 0 5 0 0 6 0 0 ...
## $ 14 : chr [1:526] "0" "0" "0" "0" ...
## $ 15 : chr [1:526] "0" "0" "0" "0" ...
## $ 16 : num [1:526] 0 0 0 0 0 20 18 90 1 0 ...
## $ 17 : num [1:526] 0 0 0 0 24 11 6 33 15 0 ...
## $ 18 : num [1:526] 0 0 0 0 1 0 1 13 0 0 ...
## $ 19 : num [1:526] 0 0 0 4 0 1 3 0 0 17 ...
## $ 20 : num [1:526] 0 0 0 0 36 0 3 0 0 21 ...
## $ 21 : num [1:526] 0 0 0 0 9 1 4 0 0 2 ...
## $ 22 : num [1:526] 0 0 62 11 0 0 5 15 0 0 ...
## $ 23 : num [1:526] 0 0 13 0 18 0 0 0 0 0 ...
## $ 24 : num [1:526] 0 0 0 0 0 0 5 0 12 0 ...
## $ 25 : num [1:526] 0 0 0 0 0 0 13 5 26 0 ...
## $ 26 : num [1:526] 0 0 0 0 145 7 2 3 2 0 ...
## $ 27 : num [1:526] 0 0 1 2 14 1 0 0 6 0 ...
## $ 28 : num [1:526] 0 0 0 6 0 1 5 1 32 0 ...
## $ 29 : num [1:526] 0 0 0 10 0 0 0 2 12 0 ...
## $ 30 : num [1:526] 0 NA 0 0 0 1 4 16 2 0 ...
## $ 31 : chr [1:526] "0" NA "0" NA ...rainfall <- rainfall %>%
mutate_if(is.character,as.numeric)## Warning: There were 3 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `14 = .Primitive("as.double")(`14`)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.rainfall## # A tibble: 526 × 33
## Year Month `1` `2` `3` `4` `5` `6` `7` `8` `9` `10` `11`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1976 1 0 0 0 0 0 0 0 0 0 0 0
## 2 1976 2 0 1 0 0 2 3 1 0 0 0 0
## 3 1976 3 0 0 0 7 0 0 0 34 0 0 0
## 4 1976 4 0 0 0 1 0 0 0 0 0 0 0
## 5 1976 5 0 0 0 0 1 12 1 24 24 76 6
## 6 1976 6 0 43 0 14 0 1 0 163 100 66 87
## 7 1976 7 0 33 10 11 4 0 0 103 0 52 28
## 8 1976 8 0 6 12 20 11 9 4 6 0 0 0
## 9 1976 9 13 0 0 0 0 0 0 0 0 3 0
## 10 1976 10 17 18 1 6 5 0 27 0 0 0 0
## # ℹ 516 more rows
## # ℹ 20 more variables: `12` <dbl>, `13` <dbl>, `14` <dbl>, `15` <dbl>,
## # `16` <dbl>, `17` <dbl>, `18` <dbl>, `19` <dbl>, `20` <dbl>, `21` <dbl>,
## # `22` <dbl>, `23` <dbl>, `24` <dbl>, `25` <dbl>, `26` <dbl>, `27` <dbl>,
## # `28` <dbl>, `29` <dbl>, `30` <dbl>, `31` <dbl>rainfall<- rainfall %>%
pivot_longer(
cols=!c( Year ,Month),
names_to = "day",
values_to = "rain"
)
rainfall## # A tibble: 16,306 × 4
## Year Month day rain
## <dbl> <dbl> <chr> <dbl>
## 1 1976 1 1 0
## 2 1976 1 2 0
## 3 1976 1 3 0
## 4 1976 1 4 0
## 5 1976 1 5 0
## 6 1976 1 6 0
## 7 1976 1 7 0
## 8 1976 1 8 0
## 9 1976 1 9 0
## 10 1976 1 10 0
## # ℹ 16,296 more rows#View(rainfall)time series
- Make a variable - average rainfall per month
library(readr)
Rainfall_dhaka <- read_csv("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/Rainfall_dhaka2.csv")## Rows: 526 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Stati, 14, 15, 31
## dbl (30): Year, Month, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 17, 18...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.library(tidyverse)
data<-Rainfall_dhaka %>%
select(!Stati) %>%
mutate_if(is.character,as.numeric) %>%
pivot_longer(cols = !c(Year,Month),names_to = "day",values_to = "rain") %>%
group_by(Year,Month) %>%
mutate(mean=mean(rain,na.rm = T))## Warning: There were 3 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `14 = .Primitive("as.double")(`14`)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.data## # A tibble: 16,306 × 5
## # Groups: Year, Month [526]
## Year Month day rain mean
## <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1976 1 1 0 0
## 2 1976 1 2 0 0
## 3 1976 1 3 0 0
## 4 1976 1 4 0 0
## 5 1976 1 5 0 0
## 6 1976 1 6 0 0
## 7 1976 1 7 0 0
## 8 1976 1 8 0 0
## 9 1976 1 9 0 0
## 10 1976 1 10 0 0
## # ℹ 16,296 more rows- Create Date variable
library(lubridate)
data<- data %>%
mutate(date=make_date(Year,Month,day))
data## # A tibble: 16,306 × 6
## # Groups: Year, Month [526]
## Year Month day rain mean date
## <dbl> <dbl> <chr> <dbl> <dbl> <date>
## 1 1976 1 1 0 0 1976-01-01
## 2 1976 1 2 0 0 1976-01-02
## 3 1976 1 3 0 0 1976-01-03
## 4 1976 1 4 0 0 1976-01-04
## 5 1976 1 5 0 0 1976-01-05
## 6 1976 1 6 0 0 1976-01-06
## 7 1976 1 7 0 0 1976-01-07
## 8 1976 1 8 0 0 1976-01-08
## 9 1976 1 9 0 0 1976-01-09
## 10 1976 1 10 0 0 1976-01-10
## # ℹ 16,296 more rows- Draw a basic line plot of rainfall data
x<-ts(data = data$mean,frequency = 365,start = c(1976 ,1,1))
plot(x)
4. Draw a line plot of average rainfall after year 2000
x4<- data%>%
filter(Year >= 2000)
plot(ts(data = x4$mean,frequency = 365,start = c(2000,1,1)))
- Decompose the time series
d<-decompose(x)- Plot the decomposed time series
plot(d)
acf(data$mean)
Econometrics 02
Lecture 05
library(AER)## Warning: package 'AER' was built under R version 4.4.2## Loading required package: car## Warning: package 'car' was built under R version 4.4.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.4.2##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 4.4.2## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.4.2##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.4.2## Loading required package: survivaldata("Journals")
#View(Journals)
Journals %>%
mutate(citrprice = price/citations)## title
## APEL Asian-Pacific Economic Literature
## SAJoEH South African Journal of Economic History
## CE Computational Economics
## MEPiTE MOCT-MOST Economic Policy in Transitional Economics
## JoSE Journal of Socio-Economics
## LabEc Labour Economics
## EDE Environment and Development Economics
## RoRPE Review of Radical Political Economics
## EoP Economics of Planning
## Mt Metroeconomica
## JoCP Journal of Consumer Policy
## REE Real Estate Economics
## DPR Development Policy Review
## MaDE Managerial and Decision Econ
## JoEF Journal of Empirical Finance
## JoEMS International Journal of Finance & Economics
## JoE&MS Journal of Economics & Management Strategy
## AEJ Atlantic Economic Journal
## EDQ Economic Development Quarterly
## CER China Economic Review
## IEP Information Economics and Policy
## AEP Australian Economic Papers
## JWE Japan and the World Economy
## JoES Journal of Economic Surveys
## IME Insurance Mathematics & Economics
## DE De Economist
## RoSE Review of Social Economy
## BP Brookings Papers
## PiRS Papers in Regional Science
## JoPE Journal of Population Economics
## FS Finance and Stochastics
## IO International Organization
## CPP Canadian Public Policy
## DvE Developing Economies
## JoHoEc Journal of Housing Economics
## JoEB Journal of Economics & Business
## JoEvEc Journal of Evolutionary Economics
## EmpEc Empirical Economics
## ER Econometric Reviews
## AgE Agricultural Economics
## QRoEF Quarterly Review of Economics & Finance
## AEL Applied Economics Letters
## EM Economic Modelling
## GC Growth and Change
## ERE Environmental and Resource Economics
## JoRE Journal of Regulatory Economics
## BoIES Bulletin of Indonesian Economic Studies
## AoRS Annals of Regional Science
## CEP Contemporary Economic Policy
## JoEcEd Journal of Economic Education
## JoJIE Journal of the Japanese & International Economies
## EaP Economics and Philosophy
## JoFI Journal of Financial Intermediation
## RaEE Resource and Energy Economics
## ORoEP Oxford Review of Economic Policy
## RoIO Review of Industrial Organization
## IRoLE International Review of Law & Economics
## JoWT Journal of World Trade
## JoM Journal of Macroeconomics
## AJoES American Journal of Economics & Sociology
## PFR Public Finance Review
## MnchSc Manchester School
## EnE Energy Economics
## ERoAE European Review of Agricultural Economics
## LH Labor History
## WA Weltwirtschaftliches Archiv / Review of World Economics
## SJoPE Scottish Journal of Political Economy
## DC Development and Change
## JoREFE Journal of Real Estate Finance and Economics
## JoPopEc Journal of Policy Modeling
## JoPKE Journal of Post-Keynesian Economics
## EP Economic Policy
## MSS Mathematical Social Sciences
## Om Omega
## HoPE History of Political Economy
## EcRec Economic Record
## JoPA Journal of Productivity Analysis
## RoIW Review of Income and Wealth
## JoRI Journal of Risk and Insurance
## TWE The World Economy
## JoAE Journal of Agricultural Economics
## EoER Economics of Education Review
## JoLR Journal of Labor Research
## TD Theory and Decision
## SCW Social Choice & Welfare
## JoCE Journal of Comparative Economics
## JoEP Journal of Economic Psychology
## IJoF International Journal of Forecasting
## Ky Kyklos
## EiEH Explorations in Economic History
## ET Economic Theory
## JoTEP Journal of Transport Economics & Policy
## BJoIR British Journal of Industrial Relations
## JoCMS Journal of Common Market Studies
## JoFM Journal of Futures Markets
## JITE Journal of Institutional & Theoretical Economics
## IJoIO International Journal of Industrial Organization
## JoFo Journal of Forecasting
## ORL Operations Research Letters
## SJoE Scandinavian Journal of Economics
## CJoE Cambridge Journal of Economics
## JoEI Journal of Economic Issues
## RSUE Regional Science & Urban Economics
## IJoGT International Journal of Game Theory
## JoRaU Journal of Risk and Uncertainty
## JoRU Journal of Financial and Quantitative Analysis
## JoDS Journal of Development Studies
## EcmtrT Econometric Theory
## JoRS Journal of Regional Science
## JoA&E Journal of Accounting & Economics
## GaEB Games and Economic Behavior
## JAE Journal of Applied Econometrics
## JoMathEc Journal of Mathematical Economics
## JoIMF Journal of International Money and Finance
## RS Regional Studies
## PDR Population & Development Review
## ES Economy and Society
## DS Decision Sciences
## JoIndEc Journal of Industrial Economics
## EE Ecological Economics
## NTJ National Tax Journal
## US Urban Studies
## HE Health Economics
## EHR Economic History Review
## RoFS Review of Financial Studies
## CJE Canadian Journal of Economics
## RA Risk Analysis
## AE Applied Economics
## OEP Oxford Economic Papers
## EDCC Economic Development & Cultural Change
## JoBF Journal of Banking and Finance
## OBoES Oxford Bulletin of Economics & Statistics
## JoEDC Journal of Economic Dynamics & Control
## SEJ Southern Economic Journal
## JoMS Journal of Management Studies
## JoBE Journal of Business Ethics
## Dm Demography
## EI Economic Inquiry
## JoEBO Journal of Economic Behavior & Organization
## JoLS Journal of Legal Studies
## JoDE Journal of Development Economics
## LanEc Land Economics
## JoLE Journal of Labor Economics
## JoEH Journal of Economic History
## JoUE Journal of Urban Economics
## Ecnmc Economica
## JoMCB Journal of Money Credit & Banking
## JoIntEc Journal of International Economics
## PC Public Choice
## RP Research Policy
## EL Economics Letters
## JoHeEc Journal of Health Economics
## JBES Journal of Business & Economic Statistics
## RJoE Rand Journal of Economics
## JoB Journal of Business
## JoHR Journal of Human Resources
## IER International Economic Review
## JoEEM Journal of Environmental Economics & Managment
## EER European Economic Review
## WD World Development
## JoPubEc Journal of Public Economics
## JEL Journal of Economic Literature
## JoLaE Journal of Law and Economics
## JEP Journal of Economic Perspectives
## AJoAE American Journal of Agricultural Economics
## JoMonEc Journal of Monetary Economics
## MngtSc Management Science
## RoE&S Review of Economics & Statistics
## RoES Review of Economic Studies
## JoE Journal of Econometrics
## JoET Journal of Economic Theory
## EJ Economic Journal
## JoFE Journal of Financial Economics
## JoCR Journal of Consumer Research
## JASA Journal of the American Statistical Association
## JoFi Journal of Finance
## QJoE Quarterly Journal of Economics
## JoPolEc Journal of Political Economy
## Ecnmt Econometrica
## AER American Economic Review
## publisher society price pages charpp citations
## APEL Blackwell no 123 440 3822 21
## SAJoEH So Afr ec history assn no 20 309 1782 22
## CE Kluwer no 443 567 2924 22
## MEPiTE Kluwer no 276 520 3234 22
## JoSE Elsevier no 295 791 3024 24
## LabEc Elsevier no 344 609 2967 24
## EDE Cambridge Univ Pres no 90 602 3185 24
## RoRPE Elsevier no 242 665 2688 27
## EoP Kluwer no 226 243 3010 28
## Mt Blackwell no 262 386 2501 30
## JoCP Kluwer no 279 578 2200 32
## REE MIT no 165 749 2496 35
## DPR Blackwell no 242 427 2731 36
## MaDE Wiley no 905 292 4472 37
## JoEF Elsevier no 355 607 3053 37
## JoEMS Springer no 375 351 4025 40
## JoE&MS MIT Press no 135 602 2394 42
## AEJ Intnl Atlantic Ec. Soc. no 171 447 3139 44
## EDQ Sage no 284 385 3318 47
## CER Elsevier no 242 167 3619 47
## IEP Elsevier no 371 442 2924 50
## AEP Blackwell no 115 495 3792 51
## JWE Elsevier no 355 577 3443 56
## JoES Blackwell no 355 674 2835 61
## IME Elsevier no 835 745 4263 61
## DE Kluwer no 223 579 2594 62
## RoSE Routledge no 172 578 2720 66
## BP Brookings Inst. no 62 394 2368 67
## PiRS Springer no 191 442 2925 68
## JoPE Springer no 411 640 3264 69
## FS Springer no 274 492 3060 69
## IO MIT no 130 814 3024 69
## CPP University of Toronto Press no 100 714 2248 75
## DvE Inst of Devel Econ no 80 539 3345 76
## JoHoEc Academic Press no 235 330 2623 78
## JoEB Elsevier no 392 563 3456 80
## JoEvEc Springer no 410 526 2989 81
## EmpEc Springer no 464 719 3366 82
## ER Dekker no 650 448 2418 85
## AgE Elsevier no 558 610 4032 87
## QRoEF Elsevier no 317 931 3036 87
## AEL Routledge no 495 774 5610 89
## EM Elsevier no 535 632 2880 89
## GC Blackwell no 123 302 2772 97
## ERE Kluwer no 717 1061 2924 97
## JoRE Kluwer no 481 628 3666 98
## BoIES ANU Press no 54 518 2858 98
## AoRS Springer no 379 580 2914 100
## CEP Oxford Univ Press no 168 566 4104 109
## JoEcEd Heldref no 82 426 2990 112
## JoJIE Academic Press no 355 451 2541 119
## EaP Cambridge Univ Pres no 95 373 2574 121
## JoFI Academic Press no 240 420 2773 121
## RaEE Elsevier no 448 416 3010 122
## ORoEP Oxford Univ Press no 255 600 4284 124
## RoIO Kluwer no 448 792 1820 126
## IRoLE Elsevier no 392 585 3479 126
## JoWT Kluwer no 475 1138 3088 134
## JoM LSU Press no 85 856 2562 137
## AJoES Blackwell no 108 578 2331 138
## PFR Sage no 394 636 2204 138
## MnchSc Blackwell no 336 618 2772 140
## EnE Elsevier no 565 495 3038 143
## ERoAE Oxford Univ Press no 255 573 3228 144
## LH Carfax no 165 515 3430 147
## WA Mohr Siebeck no 99 758 2499 150
## SJoPE Blackwell no 203 596 3036 151
## DC Blackwell no 318 888 2623 154
## JoREFE Kluwer no 476 630 3432 162
## JoPopEc Elsevier no 473 907 2160 164
## JoPKE M.E Sharpe no 186 550 2360 170
## EP Blackwell no 170 440 2856 179
## MSS Elsevier no 824 697 3420 180
## Om Elsevier no 805 780 4028 183
## HoPE Duke Univ Press no 132 767 2280 183
## EcRec Ec. Society of Australia yes 50 421 4320 188
## JoPA Kluwer no 424 567 3432 188
## RoIW Int Assn for Res in I & W yes 130 585 3312 189
## JoRI Am. Risk & Ins. Assn yes 90 720 2924 193
## TWE Blackwell no 805 1270 2842 214
## JoAE Agric. Econ. Society yes 96 277 4218 220
## EoER Elsevier no 448 471 5040 224
## JoLR George Mason Univ no 130 725 3010 227
## TD Kluwer no 595 607 2072 237
## SCW Springer no 474 674 2990 239
## JoCE Academic Press no 410 803 2442 245
## JoEP Elsevier no 395 802 2535 245
## IJoF Elsevier no 437 464 4324 249
## Ky Helbing & Lichtenhahn no 270 618 2760 254
## EiEH Academic Press no 265 447 2904 261
## ET Springer no 899 1493 2806 264
## JoTEP LSE and Univ of Bath no 133 406 3182 292
## BJoIR Blackwell no 262 663 3450 294
## JoCMS Blackwell no 506 710 2747 299
## JoFM Wiley no 1140 990 2460 302
## JITE Mohr Siebeck no 211 792 3402 313
## IJoIO Elsevier no 799 1230 3124 313
## JoFo Wiley no 760 516 4032 317
## ORL Elsevier no 442 460 4056 338
## SJoE Blackwell no 296 688 2666 339
## CJoE Oxford Univ Press no 272 814 3545 351
## JoEI Assn Ev. Economics yes 45 1003 2870 355
## RSUE Elsevier no 614 802 2967 370
## IJoGT Springer no 436 582 3366 382
## JoRaU Kluwer no 481 595 3476 383
## JoRU Univ Wash Press no 95 500 3314 388
## JoDS Frank Cass no 357 500 2604 388
## EcmtrT Cambridge Univ Pres no 280 899 2948 390
## JoRS Blackwell no 142 764 2992 397
## JoA&E Elsevier no 710 1112 2945 406
## GaEB Academic Press no 490 1197 2460 412
## JAE Wiley no 870 689 3680 412
## JoMathEc Elsevier no 1147 1340 2940 418
## JoIMF Elsevier no 743 940 2660 427
## RS Carfax no 759 911 4988 440
## PDR Population Council no 36 910 2904 445
## ES Routledge no 224 640 2970 488
## DS Georgia State Univ no 82 1067 3124 495
## JoIndEc Blackwell no 160 476 2684 497
## EE Elsevier no 1170 1990 4128 499
## NTJ National Tax Assn. no 90 785 2940 502
## US Carfax no 742 1711 3747 508
## HE Wiley no 575 736 4264 544
## EHR Blackwell no 163 846 3570 545
## RoFS Oxford Univ. Press no 175 1235 2604 547
## CJE Blackwell yes 120 1247 2992 556
## RA Kluwer no 590 1260 4628 574
## AE Routledge no 2120 1683 5445 578
## OEP Oxford Univ. Press no 205 767 2924 582
## EDCC Univ of Chicago Press no 128 889 2700 597
## JoBF Elsevier no 1539 1911 2516 602
## OBoES Blackwell no 346 545 2655 617
## JoEDC Elsevier yes 1046 1636 2945 636
## SEJ Southern Econ. Assn yes 97 1032 3680 646
## JoMS Blackwell no 686 850 3456 654
## JoBE Kluwer no 914 1270 3650 662
## Dm Pop Assn America no 85 568 6859 670
## EI Oxford Univ Press no 206 684 4134 684
## JoEBO Elsevier no 1154 1380 3330 698
## JoLS Univ of Chicago Press no 45 530 2625 700
## JoDE Elsevier no 1146 1110 2816 707
## LanEc Univ of Wisconsin Press no 95 580 3672 730
## JoLE Univ of Chicago Press no 138 600 2728 733
## JoEH Cambridge Univ Pres no 115 1200 4029 737
## JoUE Academic Press no 640 1058 2666 787
## Ecnmc Blackwell no 122 565 3120 825
## JoMCB Ohio State Univ. Press no 110 860 3168 834
## JoIntEc Elsevier no 923 1299 2898 838
## PC Kluwer no 1000 1600 2583 871
## RP Elsevier no 1234 781 4320 922
## EL Elsevier no 1492 1540 3315 930
## JoHeEc Elsevier no 810 828 2924 957
## JBES Am Stat Assn yes 90 510 6360 988
## RJoE RAND no 177 757 3450 1039
## JoB Univ of Chicago Press no 74 583 2655 1083
## JoHR Univ of Wisconsin Press no 113 837 3312 1113
## IER Blackwell no 145 969 3082 1113
## JoEEM Academic Press no 590 636 3096 1152
## EER Elsevier yes 1154 1823 2178 1243
## WD Elsevier no 1450 1145 5480 1408
## JoPubEc Elsevier no 1431 1880 2924 1437
## JEL Am Ec Assn yes 47 2632 3848 1530
## JoLaE Univ of Chicago Press no 45 850 2604 1580
## JEP Am Ec Assn yes 47 940 3036 1583
## AJoAE Am. Ag. Econ Assn no 81 1253 4368 1812
## JoMonEc Elsevier no 1010 1346 3174 1860
## MngtSc Inst for OR and MS no 334 1175 4232 2022
## RoE&S MIT press no 190 733 5600 2331
## RoES Blackwell no 180 761 3626 2411
## JoE Elsevier no 1893 1527 2178 2479
## JoET Academic Press no 1400 2000 2684 2514
## EJ Blackwell no 301 1983 3036 2540
## JoFE Elsevier no 1339 1947 2838 2676
## JoCR Univ of Chicago Press no 90 439 5336 2762
## JASA Am. Statistical Assn yes 310 1260 5664 2800
## JoFi Am. Finance Assn yes 226 2272 3036 3791
## QJoE MIT press no 148 1467 2184 4138
## JoPolEc Univ of Chicago Press no 159 1669 2640 6697
## Ecnmt Blackwell yes 178 1482 2992 7943
## AER Am Ec Assn yes 47 1867 3900 8999
## foundingyear subs field citrprice
## APEL 1986 14 General 5.857142857
## SAJoEH 1986 59 Economic History 0.909090909
## CE 1987 17 Specialized 20.136363636
## MEPiTE 1991 2 Area Studies 12.545454545
## JoSE 1972 96 Interdisciplinary 12.291666667
## LabEc 1994 15 Labor 14.333333333
## EDE 1995 14 Development 3.750000000
## RoRPE 1968 202 Specialized 8.962962963
## EoP 1987 46 Area Studies 8.071428571
## Mt 1949 46 General 8.733333333
## JoCP 1978 57 Consumer Economics 8.718750000
## REE 1973 125 Specialized 4.714285714
## DPR 1982 30 Development 6.722222222
## MaDE 1980 62 Management Science 24.459459459
## JoEF 1994 16 Finance 9.594594595
## JoEMS 1996 17 Finance 9.375000000
## JoE&MS 1992 37 Management Science 3.214285714
## AEJ 1972 148 General 3.886363636
## EDQ 1987 110 Development 6.042553191
## CER 1989 16 Area Studies 5.148936170
## IEP 1984 30 Specialized 7.420000000
## AEP 1961 61 General 2.254901961
## JWE 1988 27 Area Studies 6.339285714
## JoES 1987 45 General 5.819672131
## IME 1982 15 Insurance 13.688524590
## DE 1852 25 General 3.596774194
## RoSE 1943 203 General 2.606060606
## BP 1970 646 Public Finance 0.925373134
## PiRS 1922 59 Urban and Regional 2.808823529
## JoPE 1987 27 Demography 5.956521739
## FS 1996 31 Finance 3.971014493
## IO 1947 532 Industrial Organization 1.884057971
## CPP 1975 52 Public Finance 1.333333333
## DvE 1963 87 Development 1.052631579
## JoHoEc 1982 27 Specialized 3.012820513
## JoEB 1948 291 Business 4.900000000
## JoEvEc 1991 9 Specialized 5.061728395
## EmpEc 1976 24 General 5.658536585
## ER 1981 49 Econometrics 7.647058824
## AgE 1986 24 Agricultural Economics 6.413793103
## QRoEF 1961 323 Finance 3.643678161
## AEL 1994 69 General 5.561797753
## EM 1983 21 General 6.011235955
## GC 1969 212 Development 1.268041237
## ERE 1992 14 Natural Resources 7.391752577
## JoRE 1991 38 Industrial Organization 4.908163265
## BoIES 1964 22 Specialized 0.551020408
## AoRS 1967 77 Urban and Regional 3.790000000
## CEP 1982 290 Public Finance 1.541284404
## JoEcEd 1970 386 Specialized 0.732142857
## JoJIE 1986 50 Area Studies 2.983193277
## EaP 1984 144 Interdisciplinary 0.785123967
## JoFI 1991 48 Finance 1.983471074
## RaEE 1978 33 Natural Resources 3.672131148
## ORoEP 1985 29 Public Finance 2.056451613
## RoIO 1990 58 Industrial Organization 3.555555556
## IRoLE 1980 71 Law and Economics 3.111111111
## JoWT 1967 135 International 3.544776119
## JoM 1979 183 Macroeconomics 0.620437956
## AJoES 1941 573 Interdisciplinary 0.782608696
## PFR 1973 268 Public Finance 2.855072464
## MnchSc 1930 57 General 2.400000000
## EnE 1978 52 Natural Resources 3.951048951
## ERoAE 1974 21 Agricultural Economics 1.770833333
## LH 1960 479 Labor 1.122448980
## WA 1865 77 General 0.660000000
## SJoPE 1953 89 General 1.344370861
## DC 1970 83 Development 2.064935065
## JoREFE 1990 59 Business 2.938271605
## JoPopEc 1979 67 Specialized 2.884146341
## JoPKE 1978 225 General 1.094117647
## EP 1985 56 Public Finance 0.949720670
## MSS 1981 36 Interdisciplinary 4.577777778
## Om 1973 101 Business 4.398907104
## HoPE 1969 288 Economic History 0.721311475
## EcRec 1925 183 General 0.265957447
## JoPA 1992 25 General 2.255319149
## RoIW 1945 120 Public Finance 0.687830688
## JoRI 1964 263 Insurance 0.466321244
## TWE 1977 101 Development 3.761682243
## JoAE 1948 43 Agricultural Economics 0.436363636
## EoER 1982 73 Specialized 2.000000000
## JoLR 1980 221 Labor 0.572687225
## TD 1970 83 Theory 2.510548523
## SCW 1984 37 Public Finance 1.983263598
## JoCE 1977 185 Area Studies 1.673469388
## JoEP 1981 61 Interdisciplinary 1.612244898
## IJoF 1985 57 Specialized 1.755020080
## Ky 1948 186 General 1.062992126
## EiEH 1963 199 Economic History 1.015325670
## ET 1992 165 Theory 3.405303030
## JoTEP 1967 95 Specialized 0.455479452
## BJoIR 1963 82 Specialized 0.891156463
## JoCMS 1962 214 Area Studies 1.692307692
## JoFM 1981 152 Finance 3.774834437
## JITE 1844 33 General 0.674121406
## IJoIO 1983 81 Industrial Organization 2.552715655
## JoFo 1982 137 Specialized 2.397476341
## ORL 1982 59 Management Science 1.307692308
## SJoE 1898 98 General 0.873156342
## CJoE 1976 163 General 0.774928775
## JoEI 1967 462 Public Finance 0.126760563
## RSUE 1971 99 Urban and Regional 1.659459459
## IJoGT 1971 52 Theory 1.141361257
## JoRaU 1990 47 Insurance 1.255874674
## JoRU 1966 394 Finance 0.244845361
## JoDS 1964 227 Development 0.920103093
## EcmtrT 1985 65 Econometrics 0.717948718
## JoRS 1961 242 Urban and Regional 0.357682620
## JoA&E 1979 154 Interdisciplinary 1.748768473
## GaEB 1988 45 Theory 1.189320388
## JAE 1986 81 Econometrics 2.111650485
## JoMathEc 1974 79 Theory 2.744019139
## JoIMF 1981 120 International 1.740046838
## RS 1967 91 Urban and Regional 1.725000000
## PDR 1975 218 Demography 0.080898876
## ES 1971 97 General 0.459016393
## DS 1970 322 Management Science 0.165656566
## JoIndEc 1953 283 Industrial Organization 0.321931590
## EE 1989 40 Natural Resources 2.344689379
## NTJ 1948 552 Public Finance 0.179282869
## US 1964 222 Urban and Regional 1.460629921
## HE 1992 29 Health 1.056985294
## EHR 1947 390 Economic History 0.299082569
## RoFS 1988 95 Finance 0.319926874
## CJE 1967 305 General 0.215827338
## RA 1981 98 Insurance 1.027874564
## AE 1958 146 General 3.667820069
## OEP 1949 271 General 0.352233677
## EDCC 1952 515 Development 0.214405360
## JoBF 1977 172 Finance 2.556478405
## OBoES 1939 127 General 0.560777958
## JoEDC 1979 58 Theory 1.644654088
## SEJ 1932 524 General 0.150154799
## JoMS 1964 222 Management Science 1.048929664
## JoBE 1981 385 Business 1.380664653
## Dm 1964 413 Demography 0.126865672
## EI 1962 366 General 0.301169591
## JoEBO 1980 75 Theory 1.653295129
## JoLS 1972 238 Law and Economics 0.064285714
## JoDE 1974 142 Development 1.620933522
## LanEc 1925 437 Urban and Regional 0.130136986
## JoLE 1983 269 Labor 0.188267394
## JoEH 1941 643 Economic History 0.156037992
## JoUE 1974 230 Urban and Regional 0.813214740
## Ecnmc 1933 313 General 0.147878788
## JoMCB 1969 512 Macroeconomics 0.131894484
## JoIntEc 1971 211 International 1.101431981
## PC 1966 171 Public Finance 1.148105626
## RP 1972 34 Business 1.338394794
## EL 1978 81 General 1.604301075
## JoHeEc 1982 144 Health 0.846394984
## JBES 1983 190 Business 0.091093117
## RJoE 1970 339 Industrial Organization 0.170356112
## JoB 1928 771 Business 0.068328717
## JoHR 1966 522 Labor 0.101527403
## IER 1960 284 General 0.130278527
## JoEEM 1974 202 Natural Resources 0.512152778
## EER 1969 118 General 0.928399035
## WD 1973 160 Development 1.029829545
## JoPubEc 1972 141 Public Finance 0.995824635
## JEL 1963 972 General 0.030718954
## JoLaE 1968 542 Law and Economics 0.028481013
## JEP 1987 866 General 0.029690461
## AJoAE 1918 267 Agricultural Economics 0.044701987
## JoMonEc 1975 186 Macroeconomics 0.543010753
## MngtSc 1954 558 Management Science 0.165182987
## RoE&S 1919 523 General 0.081510082
## RoES 1933 325 General 0.074657818
## JoE 1973 129 Econometrics 0.763614361
## JoET 1969 165 Theory 0.556881464
## EJ 1890 531 General 0.118503937
## JoFE 1974 231 Finance 0.500373692
## JoCR 1974 536 Consumer Economics 0.032585083
## JASA 1971 487 Econometrics 0.110714286
## JoFi 1945 799 Finance 0.059614877
## QJoE 1886 660 General 0.035766071
## JoPolEc 1892 737 General 0.023741974
## Ecnmt 1932 346 General 0.022409669
## AER 1911 1098 General 0.005222803Journals$citrprice =Journals$price/Journals$citations
Journals$citrprice## [1] 5.857142857 0.909090909 20.136363636 12.545454545 12.291666667
## [6] 14.333333333 3.750000000 8.962962963 8.071428571 8.733333333
## [11] 8.718750000 4.714285714 6.722222222 24.459459459 9.594594595
## [16] 9.375000000 3.214285714 3.886363636 6.042553191 5.148936170
## [21] 7.420000000 2.254901961 6.339285714 5.819672131 13.688524590
## [26] 3.596774194 2.606060606 0.925373134 2.808823529 5.956521739
## [31] 3.971014493 1.884057971 1.333333333 1.052631579 3.012820513
## [36] 4.900000000 5.061728395 5.658536585 7.647058824 6.413793103
## [41] 3.643678161 5.561797753 6.011235955 1.268041237 7.391752577
## [46] 4.908163265 0.551020408 3.790000000 1.541284404 0.732142857
## [51] 2.983193277 0.785123967 1.983471074 3.672131148 2.056451613
## [56] 3.555555556 3.111111111 3.544776119 0.620437956 0.782608696
## [61] 2.855072464 2.400000000 3.951048951 1.770833333 1.122448980
## [66] 0.660000000 1.344370861 2.064935065 2.938271605 2.884146341
## [71] 1.094117647 0.949720670 4.577777778 4.398907104 0.721311475
## [76] 0.265957447 2.255319149 0.687830688 0.466321244 3.761682243
## [81] 0.436363636 2.000000000 0.572687225 2.510548523 1.983263598
## [86] 1.673469388 1.612244898 1.755020080 1.062992126 1.015325670
## [91] 3.405303030 0.455479452 0.891156463 1.692307692 3.774834437
## [96] 0.674121406 2.552715655 2.397476341 1.307692308 0.873156342
## [101] 0.774928775 0.126760563 1.659459459 1.141361257 1.255874674
## [106] 0.244845361 0.920103093 0.717948718 0.357682620 1.748768473
## [111] 1.189320388 2.111650485 2.744019139 1.740046838 1.725000000
## [116] 0.080898876 0.459016393 0.165656566 0.321931590 2.344689379
## [121] 0.179282869 1.460629921 1.056985294 0.299082569 0.319926874
## [126] 0.215827338 1.027874564 3.667820069 0.352233677 0.214405360
## [131] 2.556478405 0.560777958 1.644654088 0.150154799 1.048929664
## [136] 1.380664653 0.126865672 0.301169591 1.653295129 0.064285714
## [141] 1.620933522 0.130136986 0.188267394 0.156037992 0.813214740
## [146] 0.147878788 0.131894484 1.101431981 1.148105626 1.338394794
## [151] 1.604301075 0.846394984 0.091093117 0.170356112 0.068328717
## [156] 0.101527403 0.130278527 0.512152778 0.928399035 1.029829545
## [161] 0.995824635 0.030718954 0.028481013 0.029690461 0.044701987
## [166] 0.543010753 0.165182987 0.081510082 0.074657818 0.763614361
## [171] 0.556881464 0.118503937 0.500373692 0.032585083 0.110714286
## [176] 0.059614877 0.035766071 0.023741974 0.022409669 0.005222803plot(Journals$citrprice,Journals$subs)
plot(log(Journals$citrprice),Journals$subs)
plot(Journals$citrprice,log(Journals$subs))
plot(log(Journals$citrprice),log(Journals$subs))
# lm(y~x)
jour_lm<-lm(log(subs)~log(citrprice ),data = Journals)
jour_lm##
## Call:
## lm(formula = log(subs) ~ log(citrprice), data = Journals)
##
## Coefficients:
## (Intercept) log(citrprice)
## 4.7662 -0.5331#summary(jour_lm)
#plot(y~x)
plot(log(subs)~log(citrprice),data=Journals)
plot(log(Journals$subs)~log(Journals$citrprice))
#(plot(x,y))
plot(log(Journals$citrprice),log(Journals$subs))
abline(jour_lm)
ncitrprice<-seq(from = -6,to=4,by =0.25)
ncitrprice## [1] -6.00 -5.75 -5.50 -5.25 -5.00 -4.75 -4.50 -4.25 -4.00 -3.75 -3.50 -3.25
## [13] -3.00 -2.75 -2.50 -2.25 -2.00 -1.75 -1.50 -1.25 -1.00 -0.75 -0.50 -0.25
## [25] 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75
## [37] 3.00 3.25 3.50 3.75 4.00p<-predict( jour_lm,
newdata = data.frame(citrprice =exp(ncitrprice)),
interval = "prediction")
plot(log(Journals$citrprice),log(Journals$subs))
lines(p[,1]~ncitrprice)
lines(p[,2]~ncitrprice)
lines(p[,3]~ncitrprice)
Lecture 06
library(ggplot2)
#theme_set(theme_minimal())
# Demo dataset
head(economics)## # A tibble: 6 × 6
## date pce pop psavert uempmed unemploy
## <date> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1967-07-01 507. 198712 12.6 4.5 2944
## 2 1967-08-01 510. 198911 12.6 4.7 2945
## 3 1967-09-01 516. 199113 11.9 4.6 2958
## 4 1967-10-01 512. 199311 12.9 4.9 3143
## 5 1967-11-01 517. 199498 12.8 4.7 3066
## 6 1967-12-01 525. 199657 11.8 4.8 3018# Basic line plot
ggplot(data = economics, aes(x = date, y = pop))+
geom_line(color = "#00AFBB")
# Plot a subset of the data
ss <- subset(economics, date > as.Date("2006-1-1"))
ggplot(data = ss, aes(x = date, y = pop)) +
geom_line(color = "#FC4E07", size = 2)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#Control line size by the value of a continuous variable:
ggplot(data = economics, aes(x = date, y = pop)) +
geom_line(aes(size = unemploy/pop), color = "#FC4E07")
#Plot multiple time series data
library(tidyr)
library(dplyr)
df <- economics %>%
select(date, psavert, uempmed) %>%
gather(key = "variable", value = "value", -date)
head(df, 3)## # A tibble: 3 × 3
## date variable value
## <date> <chr> <dbl>
## 1 1967-07-01 psavert 12.6
## 2 1967-08-01 psavert 12.6
## 3 1967-09-01 psavert 11.9df2<-economics %>%
select(date, psavert, uempmed) %>%
pivot_longer(cols = !date,names_to = "variable2",values_to = "value2")
# Multiple line plot
ggplot(df, aes(x = date, y = value)) +
geom_line(aes(color = variable), size = 1) +
scale_color_manual(values = c("#00AFBB", "#E7B800")) +
theme_minimal()
ggplot(df2, aes(x = date, y = value2)) +
geom_line(aes(color = variable2), size = 1) 
# Area plot
ggplot(df, aes(x = date, y = value)) +
geom_area(aes(color = variable, fill = variable),
alpha = 0.5, position = position_dodge(0.8)) +
scale_color_manual(values = c("#00AFBB", "#E7B800")) +
scale_fill_manual(values = c("#00AFBB", "#E7B800"))
# Base plot with date axis
p <- ggplot(data = economics, aes(x = date, y = psavert)) +
geom_line(color = "#00AFBB", size = 1)
p
# Set axis limits c(min, max)
min <- as.Date("2002-1-1")
max <- NA
p + scale_x_date(limits = c(min, max))## Warning: Removed 414 rows containing missing values or values outside the scale range
## (`geom_line()`).
# Format : month/year
p + scale_x_date(date_labels = "%b/%Y")
#Add trend smoothed line
p + stat_smooth(
color = "#FC4E07", fill = "#FC4E07",
method = "loess"
)## `geom_smooth()` using formula = 'y ~ x'
## Lecture 07
library(AER)
data("Journals")
#View(Journals)
library(tidyverse)
Journals<-Journals %>%
mutate(citrprice = price/citations)
jour_lm<-lm(log(subs)~log(citrprice ),data = Journals)
jour_lm
y<-fitted(jour_lm)
jour_lm_new<-lm(log(subs)~log(citrprice )+ I(y^2)+I(y^3),data = Journals)
r2_o<-summary(jour_lm)$r.squared
r2_n<-summary(jour_lm_new)$r.squared
F_value<-((r2_n-r2_o)/2)/((1-r2_n)/(nrow(Journals)-4))
F_valuelibrary(lmtest)
#reset test
resettest(jour_lm,power = 2:3)##
## RESET test
##
## data: jour_lm
## RESET = 1.4409, df1 = 2, df2 = 176, p-value = 0.2395lm test
journals <- Journals[, c("subs", "price")]
journals$citeprice <- Journals$price/Journals$citations
jour_lm<-lm(log(subs)~log(citrprice ),data = Journals)
u <- resid(jour_lm)
x <- log(journals$citeprice)
model_u <- lm(u~x + I(x^2) + I(x^3), data = journals)
summary(model_u)##
## Call:
## lm(formula = u ~ x + I(x^2) + I(x^3), data = journals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5698 -0.4940 0.0444 0.4518 1.8269
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.067327 0.073136 0.921 0.359
## x -0.041605 0.057375 -0.725 0.469
## I(x^2) -0.024315 0.022912 -1.061 0.290
## I(x^3) 0.002377 0.007393 0.322 0.748
##
## Residual standard error: 0.7479 on 176 degrees of freedom
## Multiple R-squared: 0.01611, Adjusted R-squared: -0.0006605
## F-statistic: 0.9606 on 3 and 176 DF, p-value: 0.4127dim(journals)[1]*0.01611## [1] 2.8998F test for non-nested model
library(lmtest)
model1 <- lm(log(subs) ~ citations, data = Journals)
model2 <- lm(log(subs) ~ pages, data = Journals)
encomptest(model1,model2)## Encompassing test
##
## Model 1: log(subs) ~ citations
## Model 2: log(subs) ~ pages
## Model E: log(subs) ~ citations + pages
## Res.Df Df F Pr(>F)
## M1 vs. ME 177 -1 4.1985 0.04194 *
## M2 vs. ME 177 -1 22.5902 4.135e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Davidson Mackinnon J test
jtest(model1, model2)## J test
##
## Model 1: log(subs) ~ citations
## Model 2: log(subs) ~ pages
## Estimate Std. Error t value Pr(>|t|)
## M1 + fitted(M2) 0.44531 0.21733 2.0490 0.04194 *
## M2 + fitted(M1) 0.81201 0.17084 4.7529 4.135e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(model1)## [1] 515.3615AIC(model2)## [1] 532.7627incourse solve
AveTemp <- read_excel("~/AST 431/AveTemp.xlsx")
acf(AveTemp$AveTemp)
d<- decompose(AveTemp)
acf(diff(AveTemp$AveTemp,lag = 2))Lecture 07
#install.packages("plm")
library(plm)## Warning: package 'plm' was built under R version 4.4.2##
## Attaching package: 'plm'## The following objects are masked from 'package:dplyr':
##
## between, lag, leadlibrary(readr)
mydat <- panel_wage <- read_csv("Data/panel_wage.csv")## Rows: 4165 Columns: 22## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (22): exp, wks, occ, ind, south, smsa, ms, fem, union, ed, blk, lwage, i...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.names(mydat)## [1] "exp" "wks" "occ" "ind" "south" "smsa" "ms" "fem" "union"
## [10] "ed" "blk" "lwage" "id" "t" "tdum1" "tdum2" "tdum3" "tdum4"
## [19] "tdum5" "tdum6" "tdum7" "exp2"head(mydat)## # A tibble: 6 × 22
## exp wks occ ind south smsa ms fem union ed blk lwage id
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 32 0 0 1 0 1 0 0 9 0 5.56 1
## 2 4 43 0 0 1 0 1 0 0 9 0 5.72 1
## 3 5 40 0 0 1 0 1 0 0 9 0 6.00 1
## 4 6 39 0 0 1 0 1 0 0 9 0 6.00 1
## 5 7 42 0 1 1 0 1 0 0 9 0 6.06 1
## 6 8 35 0 1 1 0 1 0 0 9 0 6.17 1
## # ℹ 9 more variables: t <dbl>, tdum1 <dbl>, tdum2 <dbl>, tdum3 <dbl>,
## # tdum4 <dbl>, tdum5 <dbl>, tdum6 <dbl>, tdum7 <dbl>, exp2 <dbl>#Set data as panel data
pdata=pdata.frame(mydat,index=c("id","t"))
head(pdata)## exp wks occ ind south smsa ms fem union ed blk lwage id t tdum1 tdum2
## 1-1 3 32 0 0 1 0 1 0 0 9 0 5.56068 1 1 1 0
## 1-2 4 43 0 0 1 0 1 0 0 9 0 5.72031 1 2 0 1
## 1-3 5 40 0 0 1 0 1 0 0 9 0 5.99645 1 3 0 0
## 1-4 6 39 0 0 1 0 1 0 0 9 0 5.99645 1 4 0 0
## 1-5 7 42 0 1 1 0 1 0 0 9 0 6.06146 1 5 0 0
## 1-6 8 35 0 1 1 0 1 0 0 9 0 6.17379 1 6 0 0
## tdum3 tdum4 tdum5 tdum6 tdum7 exp2
## 1-1 0 0 0 0 0 9
## 1-2 0 0 0 0 0 16
## 1-3 1 0 0 0 0 25
## 1-4 0 1 0 0 0 36
## 1-5 0 0 1 0 0 49
## 1-6 0 0 0 1 0 64#Descriptive statistics
summary(pdata)## exp wks occ ind
## Min. : 1.00 Min. : 5.00 Min. :0.0000 Min. :0.0000
## 1st Qu.:11.00 1st Qu.:46.00 1st Qu.:0.0000 1st Qu.:0.0000
## Median :18.00 Median :48.00 Median :1.0000 Median :0.0000
## Mean :19.85 Mean :46.81 Mean :0.5112 Mean :0.3954
## 3rd Qu.:29.00 3rd Qu.:50.00 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :51.00 Max. :52.00 Max. :1.0000 Max. :1.0000
##
## south smsa ms fem
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:0.0000
## Median :0.0000 Median :1.0000 Median :1.0000 Median :0.0000
## Mean :0.2903 Mean :0.6538 Mean :0.8144 Mean :0.1126
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
##
## union ed blk lwage
## Min. :0.000 Min. : 4.00 Min. :0.00000 Min. :4.605
## 1st Qu.:0.000 1st Qu.:12.00 1st Qu.:0.00000 1st Qu.:6.395
## Median :0.000 Median :12.00 Median :0.00000 Median :6.685
## Mean :0.364 Mean :12.85 Mean :0.07227 Mean :6.676
## 3rd Qu.:1.000 3rd Qu.:16.00 3rd Qu.:0.00000 3rd Qu.:6.953
## Max. :1.000 Max. :17.00 Max. :1.00000 Max. :8.537
##
## id t tdum1 tdum2 tdum3
## 1 : 7 1:595 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 2 : 7 2:595 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## 3 : 7 3:595 Median :0.0000 Median :0.0000 Median :0.0000
## 4 : 7 4:595 Mean :0.1429 Mean :0.1429 Mean :0.1429
## 5 : 7 5:595 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000
## 6 : 7 6:595 Max. :1.0000 Max. :1.0000 Max. :1.0000
## (Other):4123 7:595
## tdum4 tdum5 tdum6 tdum7
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :0.0000 Median :0.0000 Median :0.0000 Median :0.0000
## Mean :0.1429 Mean :0.1429 Mean :0.1429 Mean :0.1429
## 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
##
## exp2
## Min. : 1.0
## 1st Qu.: 121.0
## Median : 324.0
## Mean : 514.4
## 3rd Qu.: 841.0
## Max. :2601.0
## #Pooled OLS estimator
pooling <- plm(lwage~exp+exp2+wks+ed , data=pdata, model="pooling")
summary(pooling)## Pooling Model
##
## Call:
## plm(formula = lwage ~ exp + exp2 + wks + ed, data = pdata, model = "pooling")
##
## Balanced Panel: n = 595, T = 7, N = 4165
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.16057670 -0.25034526 0.00027256 0.26792139 2.12969386
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 4.9080e+00 6.7330e-02 72.8945 < 2.2e-16 ***
## exp 4.4675e-02 2.3929e-03 18.6701 < 2.2e-16 ***
## exp2 -7.1563e-04 5.2794e-05 -13.5552 < 2.2e-16 ***
## wks 5.8270e-03 1.1826e-03 4.9271 8.673e-07 ***
## ed 7.6041e-02 2.2266e-03 34.1511 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 886.9
## Residual Sum of Squares: 635.41
## R-Squared: 0.28356
## Adj. R-Squared: 0.28287
## F-statistic: 411.624 on 4 and 4160 DF, p-value: < 2.22e-16#Between estimator
between <- plm(lwage~exp+exp2+wks+ed , data=pdata, model="between")
summary(between)## Oneway (individual) effect Between Model
##
## Call:
## plm(formula = lwage ~ exp + exp2 + wks + ed, data = pdata, model = "between")
##
## Balanced Panel: n = 595, T = 7, N = 4165
## Observations used in estimation: 595
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.978153 -0.220264 0.036574 0.250118 0.985629
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 4.68303917 0.21009890 22.2897 < 2.2e-16 ***
## exp 0.03815295 0.00569666 6.6974 4.953e-11 ***
## exp2 -0.00063127 0.00012568 -5.0228 6.757e-07 ***
## wks 0.01309028 0.00406592 3.2195 0.001355 **
## ed 0.07378378 0.00489848 15.0626 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 92.322
## Residual Sum of Squares: 62.187
## R-Squared: 0.32641
## Adj. R-Squared: 0.32185
## F-statistic: 71.4768 on 4 and 590 DF, p-value: < 2.22e-16# First differences estimator
firstdiff <- plm(lwage~exp+exp2+wks+ed , data=pdata, model="fd")
summary(firstdiff)## Oneway (individual) effect First-Difference Model
##
## Call:
## plm(formula = lwage ~ exp + exp2 + wks + ed, data = pdata, model = "fd")
##
## Balanced Panel: n = 595, T = 7, N = 4165
## Observations used in estimation: 3570
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.1131555 -0.0654718 -0.0095751 0.0483881 2.3295637
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 0.11706540 0.00631057 18.5507 < 2.2e-16 ***
## exp2 -0.00053212 0.00013927 -3.8207 0.0001354 ***
## wks -0.00026826 0.00056483 -0.4749 0.6348525
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 118.06
## Residual Sum of Squares: 117.58
## R-Squared: 0.004108
## Adj. R-Squared: 0.0035496
## F-statistic: 7.35691 on 2 and 3567 DF, p-value: 0.0006479# Fixed effect or within estimator
fixed <- plm(lwage~exp+exp2+wks+ed , data=pdata, model="within")
summary(fixed)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = lwage ~ exp + exp2 + wks + ed, data = pdata, model = "within")
##
## Balanced Panel: n = 595, T = 7, N = 4165
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.8120879 -0.0511128 0.0037112 0.0614250 1.9434065
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## exp 1.1379e-01 2.4689e-03 46.0888 < 2.2e-16 ***
## exp2 -4.2437e-04 5.4632e-05 -7.7678 1.036e-14 ***
## wks 8.3588e-04 5.9967e-04 1.3939 0.1634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 240.65
## Residual Sum of Squares: 82.632
## R-Squared: 0.65663
## Adj. R-Squared: 0.59916
## F-statistic: 2273.74 on 3 and 3567 DF, p-value: < 2.22e-16# Random effect estimator
random <- plm(lwage~exp+exp2+wks+ed , data=pdata, model="random")
summary(random)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = lwage ~ exp + exp2 + wks + ed, data = pdata, model = "random")
##
## Balanced Panel: n = 595, T = 7, N = 4165
##
## Effects:
## var std.dev share
## idiosyncratic 0.02317 0.15220 0.185
## individual 0.10209 0.31952 0.815
## theta: 0.8228
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.0439676 -0.1057048 0.0070992 0.1147499 2.0875839
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 3.8294e+00 9.3634e-02 40.8974 <2e-16 ***
## exp 8.8861e-02 2.8178e-03 31.5360 <2e-16 ***
## exp2 -7.7257e-04 6.2262e-05 -12.4083 <2e-16 ***
## wks 9.6577e-04 7.4329e-04 1.2993 0.1938
## ed 1.1171e-01 6.0572e-03 18.4426 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 260.94
## Residual Sum of Squares: 151.35
## R-Squared: 0.42
## Adj. R-Squared: 0.41945
## Chisq: 3012.45 on 4 DF, p-value: < 2.22e-16# LM test for random effect vs OLS
plmtest(pooling)##
## Lagrange Multiplier Test - (Honda)
##
## data: lwage ~ exp + exp2 + wks + ed
## normal = 72.056, p-value < 2.2e-16
## alternative hypothesis: significant effects# LM test(F test) for fixed effects vs OLS
pFtest(fixed, pooling)##
## F test for individual effects
##
## data: lwage ~ exp + exp2 + wks + ed
## F = 40.239, df1 = 593, df2 = 3567, p-value < 2.2e-16
## alternative hypothesis: significant effects# Hausman test for fixed versus random effects model
phtest(random, fixed)##
## Hausman Test
##
## data: lwage ~ exp + exp2 + wks + ed
## chisq = 6191.4, df = 3, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistentlibrary(foreign)
Panel <- read.dta("http://dss.princeton.edu/training/Panel101.dta")
head(Panel)## country year y y_bin x1 x2 x3 opinion op
## 1 A 1990 1342787840 1 0.2779036 -1.1079559 0.28255358 Str agree 1
## 2 A 1991 -1899660544 0 0.3206847 -0.9487200 0.49253848 Disag 0
## 3 A 1992 -11234363 0 0.3634657 -0.7894840 0.70252335 Disag 0
## 4 A 1993 2645775360 1 0.2461440 -0.8855330 -0.09439092 Disag 0
## 5 A 1994 3008334848 1 0.4246230 -0.7297683 0.94613063 Disag 0
## 6 A 1995 3229574144 1 0.4772141 -0.7232460 1.02968037 Str agree 1coplot(y ~ year|country, type = "l", data = Panel) # Lines
coplot(y ~ year|country, type = "b", data = Panel) # points and Lines
library(car)
scatterplot(y ~ year|country, boxplots =F, smooth=T, data = Panel)
scatterplot(y ~ year|country, data = Panel)
#install.packages("gplots")
library(gplots)## Warning: package 'gplots' was built under R version 4.4.2##
## Attaching package: 'gplots'## The following object is masked from 'package:stats':
##
## lowessplotmeans(y ~ country, main = "Heterogeineity across countries", data = Panel)
Regression lines for all countries
# Ensure 'country' is a factor
Panel$country <- as.factor(Panel$country)
# Compute regression lines for each country
regression_lines <- Panel %>%
group_by(country) %>%
summarize(intercept = coef(lm(y ~ x1, data = cur_data()))[1],
slope = coef(lm(y ~ x1, data = cur_data()))[2]) %>%
expand_grid(x1 = seq(min(Panel$x1), max(Panel$x1), length.out = 100)) %>%
mutate(y_pred = intercept + slope * x1)## Warning: There was 1 warning in `summarize()`.
## ℹ In argument: `intercept = coef(lm(y ~ x1, data = cur_data()))[1]`.
## ℹ In group 1: `country = A`.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.# Plot scatter points and regression lines
ggplot(data = Panel,
aes(x = x1, y = y, color = country))+
geom_point(alpha = 0.6) + # Scatter points
geom_line(data = regression_lines,
aes(x = x1, y = y_pred, color = country),
size = 1) +
labs(title = "Regression Lines for Each Country",
x = "X1",
y = "Y") +
theme_minimal()
Time Series 02
library(stats)
ac<-arima.sim(list(order=c(0,0,12),ma=c(0,rep(0,10),.9)),n=200)
acf(ac)
air<- read.csv("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/airline-passengers.csv")
ts<-ts(air,frequency = 12)
plot(ts)
ts(air)## Time Series:
## Start = 1
## End = 144
## Frequency = 1
## Month Passengers
## 1 1 112
## 2 2 118
## 3 3 132
## 4 4 129
## 5 5 121
## 6 6 135
## 7 7 148
## 8 8 148
## 9 9 136
## 10 10 119
## 11 11 104
## 12 12 118
## 13 13 115
## 14 14 126
## 15 15 141
## 16 16 135
## 17 17 125
## 18 18 149
## 19 19 170
## 20 20 170
## 21 21 158
## 22 22 133
## 23 23 114
## 24 24 140
## 25 25 145
## 26 26 150
## 27 27 178
## 28 28 163
## 29 29 172
## 30 30 178
## 31 31 199
## 32 32 199
## 33 33 184
## 34 34 162
## 35 35 146
## 36 36 166
## 37 37 171
## 38 38 180
## 39 39 193
## 40 40 181
## 41 41 183
## 42 42 218
## 43 43 230
## 44 44 242
## 45 45 209
## 46 46 191
## 47 47 172
## 48 48 194
## 49 49 196
## 50 50 196
## 51 51 236
## 52 52 235
## 53 53 229
## 54 54 243
## 55 55 264
## 56 56 272
## 57 57 237
## 58 58 211
## 59 59 180
## 60 60 201
## 61 61 204
## 62 62 188
## 63 63 235
## 64 64 227
## 65 65 234
## 66 66 264
## 67 67 302
## 68 68 293
## 69 69 259
## 70 70 229
## 71 71 203
## 72 72 229
## 73 73 242
## 74 74 233
## 75 75 267
## 76 76 269
## 77 77 270
## 78 78 315
## 79 79 364
## 80 80 347
## 81 81 312
## 82 82 274
## 83 83 237
## 84 84 278
## 85 85 284
## 86 86 277
## 87 87 317
## 88 88 313
## 89 89 318
## 90 90 374
## 91 91 413
## 92 92 405
## 93 93 355
## 94 94 306
## 95 95 271
## 96 96 306
## 97 97 315
## 98 98 301
## 99 99 356
## 100 100 348
## 101 101 355
## 102 102 422
## 103 103 465
## 104 104 467
## 105 105 404
## 106 106 347
## 107 107 305
## 108 108 336
## 109 109 340
## 110 110 318
## 111 111 362
## 112 112 348
## 113 113 363
## 114 114 435
## 115 115 491
## 116 116 505
## 117 117 404
## 118 118 359
## 119 119 310
## 120 120 337
## 121 121 360
## 122 122 342
## 123 123 406
## 124 124 396
## 125 125 420
## 126 126 472
## 127 127 548
## 128 128 559
## 129 129 463
## 130 130 407
## 131 131 362
## 132 132 405
## 133 133 417
## 134 134 391
## 135 135 419
## 136 136 461
## 137 137 472
## 138 138 535
## 139 139 622
## 140 140 606
## 141 141 508
## 142 142 461
## 143 143 390
## 144 144 432decompose(ts)## $x
## Month Passengers
## Jan 1 1 112
## Feb 1 2 118
## Mar 1 3 132
## Apr 1 4 129
## May 1 5 121
## Jun 1 6 135
## Jul 1 7 148
## Aug 1 8 148
## Sep 1 9 136
## Oct 1 10 119
## Nov 1 11 104
## Dec 1 12 118
## Jan 2 13 115
## Feb 2 14 126
## Mar 2 15 141
## Apr 2 16 135
## May 2 17 125
## Jun 2 18 149
## Jul 2 19 170
## Aug 2 20 170
## Sep 2 21 158
## Oct 2 22 133
## Nov 2 23 114
## Dec 2 24 140
## Jan 3 25 145
## Feb 3 26 150
## Mar 3 27 178
## Apr 3 28 163
## May 3 29 172
## Jun 3 30 178
## Jul 3 31 199
## Aug 3 32 199
## Sep 3 33 184
## Oct 3 34 162
## Nov 3 35 146
## Dec 3 36 166
## Jan 4 37 171
## Feb 4 38 180
## Mar 4 39 193
## Apr 4 40 181
## May 4 41 183
## Jun 4 42 218
## Jul 4 43 230
## Aug 4 44 242
## Sep 4 45 209
## Oct 4 46 191
## Nov 4 47 172
## Dec 4 48 194
## Jan 5 49 196
## Feb 5 50 196
## Mar 5 51 236
## Apr 5 52 235
## May 5 53 229
## Jun 5 54 243
## Jul 5 55 264
## Aug 5 56 272
## Sep 5 57 237
## Oct 5 58 211
## Nov 5 59 180
## Dec 5 60 201
## Jan 6 61 204
## Feb 6 62 188
## Mar 6 63 235
## Apr 6 64 227
## May 6 65 234
## Jun 6 66 264
## Jul 6 67 302
## Aug 6 68 293
## Sep 6 69 259
## Oct 6 70 229
## Nov 6 71 203
## Dec 6 72 229
## Jan 7 73 242
## Feb 7 74 233
## Mar 7 75 267
## Apr 7 76 269
## May 7 77 270
## Jun 7 78 315
## Jul 7 79 364
## Aug 7 80 347
## Sep 7 81 312
## Oct 7 82 274
## Nov 7 83 237
## Dec 7 84 278
## Jan 8 85 284
## Feb 8 86 277
## Mar 8 87 317
## Apr 8 88 313
## May 8 89 318
## Jun 8 90 374
## Jul 8 91 413
## Aug 8 92 405
## Sep 8 93 355
## Oct 8 94 306
## Nov 8 95 271
## Dec 8 96 306
## Jan 9 97 315
## Feb 9 98 301
## Mar 9 99 356
## Apr 9 100 348
## May 9 101 355
## Jun 9 102 422
## Jul 9 103 465
## Aug 9 104 467
## Sep 9 105 404
## Oct 9 106 347
## Nov 9 107 305
## Dec 9 108 336
## Jan 10 109 340
## Feb 10 110 318
## Mar 10 111 362
## Apr 10 112 348
## May 10 113 363
## Jun 10 114 435
## Jul 10 115 491
## Aug 10 116 505
## Sep 10 117 404
## Oct 10 118 359
## Nov 10 119 310
## Dec 10 120 337
## Jan 11 121 360
## Feb 11 122 342
## Mar 11 123 406
## Apr 11 124 396
## May 11 125 420
## Jun 11 126 472
## Jul 11 127 548
## Aug 11 128 559
## Sep 11 129 463
## Oct 11 130 407
## Nov 11 131 362
## Dec 11 132 405
## Jan 12 133 417
## Feb 12 134 391
## Mar 12 135 419
## Apr 12 136 461
## May 12 137 472
## Jun 12 138 535
## Jul 12 139 622
## Aug 12 140 606
## Sep 12 141 508
## Oct 12 142 461
## Nov 12 143 390
## Dec 12 144 432
##
## $seasonal
## Jan Feb Mar Apr May Jun Jul
## 1 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 2 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 3 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 4 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 5 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 6 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 7 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 8 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 9 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 10 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 11 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 12 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 13 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 14 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 15 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 16 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 17 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 18 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 19 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 20 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 21 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 22 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 23 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## 24 -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389 31.915404
## Aug Sep Oct Nov Dec
## 1 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 2 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 3 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 4 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 5 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 6 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 7 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 8 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 9 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 10 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 11 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 12 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 13 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 14 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 15 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 16 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 17 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 18 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 19 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 20 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 21 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 22 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 23 31.411616 8.260101 -10.321338 -26.796717 -14.309975
## 24 31.411616 8.260101 -10.321338 -26.796717 -14.309975
##
## $trend
## [,1] [,2]
## Jan 1 NA NA
## Feb 1 NA NA
## Mar 1 NA NA
## Apr 1 NA NA
## May 1 NA NA
## Jun 1 NA NA
## Jul 1 7 126.7917
## Aug 1 8 127.2500
## Sep 1 9 127.9583
## Oct 1 10 128.5833
## Nov 1 11 129.0000
## Dec 1 12 129.7500
## Jan 2 13 131.2500
## Feb 2 14 133.0833
## Mar 2 15 134.9167
## Apr 2 16 136.4167
## May 2 17 137.4167
## Jun 2 18 138.7500
## Jul 2 19 140.9167
## Aug 2 20 143.1667
## Sep 2 21 145.7083
## Oct 2 22 148.4167
## Nov 2 23 151.5417
## Dec 2 24 154.7083
## Jan 3 25 157.1250
## Feb 3 26 159.5417
## Mar 3 27 161.8333
## Apr 3 28 164.1250
## May 3 29 166.6667
## Jun 3 30 169.0833
## Jul 3 31 171.2500
## Aug 3 32 173.5833
## Sep 3 33 175.4583
## Oct 3 34 176.8333
## Nov 3 35 178.0417
## Dec 3 36 180.1667
## Jan 4 37 183.1250
## Feb 4 38 186.2083
## Mar 4 39 189.0417
## Apr 4 40 191.2917
## May 4 41 193.5833
## Jun 4 42 195.8333
## Jul 4 43 198.0417
## Aug 4 44 199.7500
## Sep 4 45 202.2083
## Oct 4 46 206.2500
## Nov 4 47 210.4167
## Dec 4 48 213.3750
## Jan 5 49 215.8333
## Feb 5 50 218.5000
## Mar 5 51 220.9167
## Apr 5 52 222.9167
## May 5 53 224.0833
## Jun 5 54 224.7083
## Jul 5 55 225.3333
## Aug 5 56 225.3333
## Sep 5 57 224.9583
## Oct 5 58 224.5833
## Nov 5 59 224.4583
## Dec 5 60 225.5417
## Jan 6 61 228.0000
## Feb 6 62 230.4583
## Mar 6 63 232.2500
## Apr 6 64 233.9167
## May 6 65 235.6250
## Jun 6 66 237.7500
## Jul 6 67 240.5000
## Aug 6 68 243.9583
## Sep 6 69 247.1667
## Oct 6 70 250.2500
## Nov 6 71 253.5000
## Dec 6 72 257.1250
## Jan 7 73 261.8333
## Feb 7 74 266.6667
## Mar 7 75 271.1250
## Apr 7 76 275.2083
## May 7 77 278.5000
## Jun 7 78 281.9583
## Jul 7 79 285.7500
## Aug 7 80 289.3333
## Sep 7 81 293.2500
## Oct 7 82 297.1667
## Nov 7 83 301.0000
## Dec 7 84 305.4583
## Jan 8 85 309.9583
## Feb 8 86 314.4167
## Mar 8 87 318.6250
## Apr 8 88 321.7500
## May 8 89 324.5000
## Jun 8 90 327.0833
## Jul 8 91 329.5417
## Aug 8 92 331.8333
## Sep 8 93 334.4583
## Oct 8 94 337.5417
## Nov 8 95 340.5417
## Dec 8 96 344.0833
## Jan 9 97 348.2500
## Feb 9 98 353.0000
## Mar 9 99 357.6250
## Apr 9 100 361.3750
## May 9 101 364.5000
## Jun 9 102 367.1667
## Jul 9 103 369.4583
## Aug 9 104 371.2083
## Sep 9 105 372.1667
## Oct 9 106 372.4167
## Nov 9 107 372.7500
## Dec 9 108 373.6250
## Jan 10 109 375.2500
## Feb 10 110 377.9167
## Mar 10 111 379.5000
## Apr 10 112 380.0000
## May 10 113 380.7083
## Jun 10 114 380.9583
## Jul 10 115 381.8333
## Aug 10 116 383.6667
## Sep 10 117 386.5000
## Oct 10 118 390.3333
## Nov 10 119 394.7083
## Dec 10 120 398.6250
## Jan 11 121 402.5417
## Feb 11 122 407.1667
## Mar 11 123 411.8750
## Apr 11 124 416.3333
## May 11 125 420.5000
## Jun 11 126 425.5000
## Jul 11 127 430.7083
## Aug 11 128 435.1250
## Sep 11 129 437.7083
## Oct 11 130 440.9583
## Nov 11 131 445.8333
## Dec 11 132 450.6250
## Jan 12 133 456.3333
## Feb 12 134 461.3750
## Mar 12 135 465.2083
## Apr 12 136 469.3333
## May 12 137 472.7500
## Jun 12 138 475.0417
## Jul 12 NA NA
## Aug 12 NA NA
## Sep 12 NA NA
## Oct 12 NA NA
## Nov 12 NA NA
## Dec 12 NA NA
##
## $random
## x - seasonal.x.Month x - seasonal.x.Passengers
## Jan 1 NA NA
## Feb 1 NA NA
## Mar 1 NA NA
## Apr 1 NA NA
## May 1 NA NA
## Jun 1 NA NA
## Jul 1 -31.915404 -10.70707071
## Aug 1 -31.411616 -10.66161616
## Sep 1 -8.260101 -0.21843434
## Oct 1 10.321338 0.73800505
## Nov 1 26.796717 1.79671717
## Dec 1 14.309975 2.55997475
## Jan 2 12.374369 -3.87563131
## Feb 2 18.094066 11.01073232
## Mar 2 1.120581 7.20391414
## Apr 2 4.018308 2.60164141
## May 2 2.253157 -10.16351010
## Jun 2 -17.701389 -7.45138889
## Jul 2 -31.915404 -2.83207071
## Aug 2 -31.411616 -4.57828283
## Sep 2 -8.260101 4.03156566
## Oct 2 10.321338 -5.09532828
## Nov 2 26.796717 -10.74494949
## Dec 2 14.309975 -0.39835859
## Jan 3 12.374369 0.24936869
## Feb 3 18.094066 8.55239899
## Mar 3 1.120581 17.28724747
## Apr 3 4.018308 2.89330808
## May 3 2.253157 7.58648990
## Jun 3 -17.701389 -8.78472222
## Jul 3 -31.915404 -4.16540404
## Aug 3 -31.411616 -5.99494949
## Sep 3 -8.260101 0.28156566
## Oct 3 10.321338 -4.51199495
## Nov 3 26.796717 -5.24494949
## Dec 3 14.309975 0.14330808
## Jan 4 12.374369 0.24936869
## Feb 4 18.094066 11.88573232
## Mar 4 1.120581 5.07891414
## Apr 4 4.018308 -6.27335859
## May 4 2.253157 -8.33017677
## Jun 4 -17.701389 4.46527778
## Jul 4 -31.915404 0.04292929
## Aug 4 -31.411616 10.83838384
## Sep 4 -8.260101 -1.46843434
## Oct 4 10.321338 -4.92866162
## Nov 4 26.796717 -11.61994949
## Dec 4 14.309975 -5.06502525
## Jan 5 12.374369 -7.45896465
## Feb 5 18.094066 -4.40593434
## Mar 5 1.120581 16.20391414
## Apr 5 4.018308 16.10164141
## May 5 2.253157 7.16982323
## Jun 5 -17.701389 0.59027778
## Jul 5 -31.915404 6.75126263
## Aug 5 -31.411616 15.25505051
## Sep 5 -8.260101 3.78156566
## Oct 5 10.321338 -3.26199495
## Nov 5 26.796717 -17.66161616
## Dec 5 14.309975 -10.23169192
## Jan 6 12.374369 -11.62563131
## Feb 6 18.094066 -24.36426768
## Mar 6 1.120581 3.87058081
## Apr 6 4.018308 -2.89835859
## May 6 2.253157 0.62815657
## Jun 6 -17.701389 8.54861111
## Jul 6 -31.915404 29.58459596
## Aug 6 -31.411616 17.63005051
## Sep 6 -8.260101 3.57323232
## Oct 6 10.321338 -10.92866162
## Nov 6 26.796717 -23.70328283
## Dec 6 14.309975 -13.81502525
## Jan 7 12.374369 -7.45896465
## Feb 7 18.094066 -15.57260101
## Mar 7 1.120581 -3.00441919
## Apr 7 4.018308 -2.19002525
## May 7 2.253157 -6.24684343
## Jun 7 -17.701389 15.34027778
## Jul 7 -31.915404 46.33459596
## Aug 7 -31.411616 26.25505051
## Sep 7 -8.260101 10.48989899
## Oct 7 10.321338 -12.84532828
## Nov 7 26.796717 -37.20328283
## Dec 7 14.309975 -13.14835859
## Jan 8 12.374369 -13.58396465
## Feb 8 18.094066 -19.32260101
## Mar 8 1.120581 -0.50441919
## Apr 8 4.018308 -4.73169192
## May 8 2.253157 -4.24684343
## Jun 8 -17.701389 29.21527778
## Jul 8 -31.915404 51.54292929
## Aug 8 -31.411616 41.75505051
## Sep 8 -8.260101 12.28156566
## Oct 8 10.321338 -21.22032828
## Nov 8 26.796717 -42.74494949
## Dec 8 14.309975 -23.77335859
## Jan 9 12.374369 -20.87563131
## Feb 9 18.094066 -33.90593434
## Mar 9 1.120581 -0.50441919
## Apr 9 4.018308 -9.35669192
## May 9 2.253157 -7.24684343
## Jun 9 -17.701389 37.13194444
## Jul 9 -31.915404 63.62626263
## Aug 9 -31.411616 64.38005051
## Sep 9 -8.260101 23.57323232
## Oct 9 10.321338 -15.09532828
## Nov 9 26.796717 -40.95328283
## Dec 9 14.309975 -23.31502525
## Jan 10 12.374369 -22.87563131
## Feb 10 18.094066 -41.82260101
## Mar 10 1.120581 -16.37941919
## Apr 10 4.018308 -27.98169192
## May 10 2.253157 -15.45517677
## Jun 10 -17.701389 36.34027778
## Jul 10 -31.915404 77.25126263
## Aug 10 -31.411616 89.92171717
## Sep 10 -8.260101 9.23989899
## Oct 10 10.321338 -21.01199495
## Nov 10 26.796717 -57.91161616
## Dec 10 14.309975 -47.31502525
## Jan 11 12.374369 -30.16729798
## Feb 11 18.094066 -47.07260101
## Mar 11 1.120581 -4.75441919
## Apr 11 4.018308 -16.31502525
## May 11 2.253157 1.75315657
## Jun 11 -17.701389 28.79861111
## Jul 11 -31.915404 85.37626263
## Aug 11 -31.411616 92.46338384
## Sep 11 -8.260101 17.03156566
## Oct 11 10.321338 -23.63699495
## Nov 11 26.796717 -57.03661616
## Dec 11 14.309975 -31.31502525
## Jan 12 12.374369 -26.95896465
## Feb 12 18.094066 -52.28093434
## Mar 12 1.120581 -45.08775253
## Apr 12 4.018308 -4.31502525
## May 12 2.253157 1.50315657
## Jun 12 -17.701389 42.25694444
## Jul 12 NA NA
## Aug 12 NA NA
## Sep 12 NA NA
## Oct 12 NA NA
## Nov 12 NA NA
## Dec 12 NA NA
##
## $figure
## [1] -12.374369 -18.094066 -1.120581 -4.018308 -2.253157 17.701389
## [7] 31.915404 31.411616 8.260101 -10.321338 -26.796717 -14.309975
##
## $type
## [1] "additive"
##
## attr(,"class")
## [1] "decomposed.ts"MA(1) process positive coefficient
ma.sim1 <- arima.sim(list(order = c(0,0,1), ma = 0.9), n = 100)
plot(ma.sim1,
ylab=expression(Y[t]),
xlab="Time",
type="o",
main="MA(1) simulation")
plot(ma.sim1,type="o")
acf(ma.sim1,main="Sample ACF")
MA(1) process negative coefficient
ma.sim2 <- arima.sim(list(order = c(0,0,1), ma = -0.9), n = 100)
plot(ma.sim2,
ylab=expression(Y[t]),
xlab="Time",
type="o",
main="MA(1) simulation")
acf(ma.sim2,main="Sample ACF")
MA(2) process
ma.sim3 <- arima.sim(list(order = c(0,0,2), ma = c(-0.9,0.7)), n = 100)
plot(ma.sim3,
ylab=expression(Y[t]),
xlab="Time",
type="o",
main="MA(2) simulation")
acf(ma.sim3,main="Sample ACF")[1:20] # gives values
##
## Autocorrelations of series 'ma.sim3', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## -0.734 0.401 -0.140 0.078 -0.102 0.079 -0.129 0.141 -0.144 0.128 -0.076
## 12 13 14 15 16 17 18 19 20
## 0.037 0.031 -0.060 0.044 -0.035 0.000 -0.055 0.094 -0.126plot(acf(ma.sim3, plot = F)[1:20]) # plots
AR(1) process positive coefficient
ar.sim.1 <- arima.sim(list(order = c(1,0,0), ar = 0.9), n = 100)
plot(ar.sim.1,
ylab=expression(Y[t]),
xlab="Time",
type="o",
main="AR(1) simulation")
acf(ar.sim.1,main="Sample ACF")
pacf(ar.sim.1)
# AR(1) process negative coefficient
ar.sim.2 <- arima.sim(list(order = c(1,0,0), ar = -0.9), n = 100)
plot(ar.sim.2,
ylab=expression(Y[t]),
xlab="Time",
type="o",
main="AR(1) simulation")
acf(ar.sim.2,main="Sample ACF")
pacf(ar.sim.2)
#AR(2) process
ar.sim.3 <- arima.sim(list(order = c(2,0,0), ar = c(0.7,-0.9)), n = 100)
plot(ar.sim.3,
ylab=expression(Y[t]),
xlab="Time",
type="o",
main="AR(2) simulation")
acf(ar.sim.3,main="Sample ACF")[1:20]
##
## Autocorrelations of series 'ar.sim.3', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.347 -0.656 -0.729 0.103 0.689 0.369 -0.362 -0.565 -0.049 0.454 0.316
## 12 13 14 15 16 17 18 19 20
## -0.197 -0.382 -0.071 0.257 0.196 -0.085 -0.187 -0.033 0.091pacf(ar.sim.3)
#ARMA(1,1) process
X<-arima.sim(n=1000,list(ar=c(0.5),ma=c(0.2)))
acf(X)
pacf(X)
# ARMA(1,1) process (different coefficient)
Y<-arima.sim(n=1000,list(order=c(1,0,1),ar=0.9,ma=0.2))
acf(Y)
pacf(Y)
Fitting AR(1) on raw data
data <- c(14.3, 14.6, 13.5, 14.2, 12.1, 14.19, 14.55, 13.6, 14.59, 16.6, 15.4, 12.89, 12.2, 12.15, 10.89, 11.11, 11.98, 13.44, 14.05, 13.91, 14.1, 12.66, 14.6, 16.7, 15.4, 17.1, 15.45, 16.76, 15.7, 14.1, 15.3, 16.4, 17.09, 16.8, 16.98, 16.39, 15.8, 15.1, 13.7, 13.79, 15.7)
acf(data)
pacf(data)
a <- arima(data, c(1,0,0))
names(a)## [1] "coef" "sigma2" "var.coef" "mask" "loglik" "aic"
## [7] "arma" "residuals" "call" "series" "code" "n.cond"
## [13] "nobs" "model"acf(a$residuals)
pacf(a$residuals)
Using co2_mm_mlo.csv data
Generate a white noise process
WN <- arima.sim(model = list(order = c(0, 0, 0)), n = 200)
plot(WN, main="White Noise Series")
library(readr)
co2_csv <- read_csv("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/co2_mm_mlo.csv")## Rows: 761 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (8): year, month, decimal date, average, deseasonalized, ndays, sdev, unc
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#View(co2_mm_mlo)
co2 <- co2_csv$average
#date_co2 <- timeSequence(from="1958-03-01",to="2021-07-30", by="month")
#co2.ts <- timeSeries(data=co2,date_co2, units="co2")
#plot(co2.ts, main = "Time Series Plot of co2")
#my code
co2<-ts(co2,start = c(1958,3),frequency = 12)
plot(co2)
Plotting various components of the Time Series
co2.stl<-stl(co2,
s.window =12,
t.window=0.75*length(co2),
t.degree=0)
plot(co2.stl$time.series[,"trend"], main = "Trend")
plot(co2.stl$time.series[,"seasonal"], main = "Seasonal")
plot(co2.stl$time.series[,"remainder"], main = "Remainder")
co2d <- diff(co2.stl$time.series[,"remainder"])
plot(co2d)
#Checking PACF and ACF
pacf(as.vector(co2d),main="PACF Plot") 
acf(co2d, main="ACF Plot")
#Fitting an autoregressive model
co2dar<-ar(co2d)
str(co2dar)## List of 15
## $ order : int 15
## $ ar : num [1:15] -0.2999 -0.0914 -0.0512 -0.0185 0.0281 ...
## $ var.pred : num 0.0952
## $ x.mean : num 0.0631
## $ aic : Named num [1:29] 56.77 8.97 9.4 11.33 13.32 ...
## ..- attr(*, "names")= chr [1:29] "0" "1" "2" "3" ...
## $ n.used : int 760
## $ n.obs : int 760
## $ order.max : num 28
## $ partialacf : num [1:28, 1, 1] -0.25183 -0.04548 -0.00961 0.00431 0.0396 ...
## $ resid : Time-Series [1:760] from 1958 to 2022: NA NA NA NA NA NA NA NA NA NA ...
## $ method : chr "Yule-Walker"
## $ series : chr "co2d"
## $ frequency : num 12
## $ call : language ar(x = co2d)
## $ asy.var.coef: num [1:15, 1:15] 1.34e-03 3.94e-04 1.15e-04 5.63e-05 1.43e-05 ...
## - attr(*, "class")= chr "ar"summary(co2dar)## Length Class Mode
## order 1 -none- numeric
## ar 15 -none- numeric
## var.pred 1 -none- numeric
## x.mean 1 -none- numeric
## aic 29 -none- numeric
## n.used 1 -none- numeric
## n.obs 1 -none- numeric
## order.max 1 -none- numeric
## partialacf 28 -none- numeric
## resid 760 ts numeric
## method 1 -none- character
## series 1 -none- character
## frequency 1 -none- numeric
## call 2 -none- call
## asy.var.coef 225 -none- numeric#Fitting ARIMA(1,0,1) to this data
fit<-arima(co2d,order=c(1,0,1))
tsdiag(fit)
qqnorm(fit$residuals)
# 1.1 Reading and Cleaning the Rainfall in Dhaka data
library(readr)
Rainfall_dhaka2 <- read_csv("Data/Rainfall_dhaka2.csv")## Rows: 526 Columns: 34
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Stati, 14, 15, 31
## dbl (30): Year, Month, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16, 17, 18...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rainfall<-Rainfall_dhaka2 %>%
mutate_at(vars("14","15","31"),as.integer)## Warning: There were 3 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `14 = .Primitive("as.integer")(`14`)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2 remaining warnings.# 1.2 Makingavariable:Averagerainfallpermonth
Rainfall$avg<-rowMeans(Rainfall[,4:ncol(Rainfall)],na.rm = T)
Rain<-Rainfall[,c(2,3,35)]
# 1.3 Abasic line plot of average rainfall per month (1976-2019)
# Preparing the Time Series Data
Rain_ts<-ts(data = Rain$avg,start = c(1976,1),frequency = 12)
plot(Rain_ts)
# 1.4 Checking out the individual components
rain.stl<-stl(Rain_ts,s.window = 12,t.window = 0.75*length(Rain_ts),t.degree = 0)
#Trend
plot(rain.stl$time.series[,"trend"])
#Seasonality
plot(rain.stl$time.series[,"seasonal"])
# Remainder
plot(rain.stl$time.series[,"remainder"])
# 1.5 Checking ACF
#Remainder
rain_res<-rain.stl$time.series[,"remainder"]
acf(rain_res)
# 1.6 Checking PACF
pacf(rain_res)
# 1.7 Fitting Models
# Wefit some closely related models that might best describe this time series.
m1 <- arima(rain_res,order = c(0,0,1))
m2 <- arima(rain_res,order = c(1,0,0))
m3 <- arima(rain_res,order = c(1,1,1))
m4 <- arima(rain_res,order = c(0,0,2))
m5 <- arima(rain_res,order = c(1,0,2))
m6 <- arima(rain_res,order = c(1,1,2))
m7 <- arima(rain_res,order = c(2,0,2))
m8 <- arima(rain_res,order = c(2,1,2))
m9 <- arima(rain_res,order = c(0,0,3))
# 1.8 Comparing AIC values for different models
AIC(m1)## [1] 2689.181 AIC(m2)## [1] 2686.034 AIC(m3)## [1] 2687.775 AIC(m4)## [1] 2684.678 AIC(m5)## [1] 2686.675 AIC(m6)## [1] 2687.325 AIC(m7)## [1] 2688.624 AIC(m8)## [1] 2688.275 AIC(m9)## [1] 2686.675 # 1.9 Summary estimates of ARIMA(0,0,2)
broom::tidy(m4)## # A tibble: 3 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 ma1 0.180 0.0433
## 2 ma2 0.110 0.0428
## 3 intercept -0.0240 0.173library(zoo)
library(forecast)## Warning: package 'forecast' was built under R version 4.4.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(Metrics)## Warning: package 'Metrics' was built under R version 4.4.3##
## Attaching package: 'Metrics'## The following object is masked from 'package:forecast':
##
## accuracy#Seasonal MA(12) process
arima.m<-arima.sim(list(order = c(0,0,12), ma = c(0.7,rep(0,10),0.9)), n = 200)
acf(arima.m)
# A custom seasonal AR(13) process
ar_pacf<-ARMAacf (ar = c(.6,0,0,0,0,0,0,0,0,0,0,.5,-.30),lag.max=30,pacf=T)
plot(ar_pacf, type='h')
library(readr)
airline_passengers <- read_csv("Data/airline-passengers.csv")## Rows: 144 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Month
## dbl (1): Passengers
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#View(airline_passengers)
airline_passengers$Month <- as.yearmon(airline_passengers$Month)
max(airline_passengers$Passengers)## [1] 622airline_passengers[which.max(airline_passengers$Passengers), ]## # A tibble: 1 × 2
## Month Passengers
## <yearmon> <dbl>
## 1 Jul 1960 622Divides the dataframe into train and test sets and then tries to predict the test set using the trained data. It is a naive approach. Calculates monthly averages as future predictions
test_set_start_date <- as.yearmon("Jan 1959")
train_set <- subset(airline_passengers, Month < test_set_start_date)
test_set <- subset(airline_passengers, Month >= test_set_start_date)
train_set <- subset(airline_passengers, Month < "1959-01")
test_set <- subset(airline_passengers, Month >= "1959-01")
dim(train_set)## [1] 120 2dim(test_set)## [1] 24 2train_set_avg <- mean(train_set$Passengers)
simple_avg_predictions <- data.frame(
Month = test_set$Month,
Passengers = rep(train_set_avg, nrow(test_set))
)
ggplot() +
geom_line(data = train_set, aes(x = Month, y = Passengers, color = "Training"), size = 1) +
geom_line(data = test_set, aes(x = Month, y = Passengers, color = "Testing"), size = 1) +
labs(
title = "Airline Passengers - Training and Testing Sets",
x = "Date",
y = "Passengers"
) +
scale_color_manual(values = c("Training" = "#12355B", "Testing" = "#D72638"), name = "Airline passengers") +
theme_minimal() +
theme(plot.title = element_text(size = 20))## Warning: The `trans` argument of `continuous_scale()` is deprecated as of ggplot2 3.5.0.
## ℹ Please use the `transform` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplot()+
geom_line(data = train_set,aes(x = Month, y = Passengers,colour = "training"))+
geom_line(data = test_set,aes(x = Month, y = Passengers,colour = "testing"))
Converts the training set to a time series data, then fits a seasonal ARIMA(1,0,0)(2,1,0)[12].
series_ts <- ts(train_set$Passengers, frequency = 12)
series_ts## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1 112 118 132 129 121 135 148 148 136 119 104 118
## 2 115 126 141 135 125 149 170 170 158 133 114 140
## 3 145 150 178 163 172 178 199 199 184 162 146 166
## 4 171 180 193 181 183 218 230 242 209 191 172 194
## 5 196 196 236 235 229 243 264 272 237 211 180 201
## 6 204 188 235 227 234 264 302 293 259 229 203 229
## 7 242 233 267 269 270 315 364 347 312 274 237 278
## 8 284 277 317 313 318 374 413 405 355 306 271 306
## 9 315 301 356 348 355 422 465 467 404 347 305 336
## 10 340 318 362 348 363 435 491 505 404 359 310 337mod2<-Arima(series_ts,order=c(1, 0, 0),
seasonal=list(order=c(2, 1, 0), period=12))
mod2## Series: series_ts
## ARIMA(1,0,0)(2,1,0)[12]
##
## Coefficients:
## ar1 sar1 sar2
## 0.9418 -0.0637 0.0781
## s.e. 0.0313 0.1079 0.1107
##
## sigma^2 = 107.1: log likelihood = -405.29
## AIC=818.58 AICc=818.97 BIC=829.31Seasonal ARIMA (Cont…)
This lecture is based in sarima(1).
Generates a hypothetical sales data that has trend, seasonality, and random white noise, based on a monthly data starting from 2020 Jan 1st.
library(forecast)
library(tseries)## Warning: package 'tseries' was built under R version 4.4.3set.seed(123) # For reproducibility
months <- seq(as.Date("2020-01-01"),
by = "month",
length.out = 36)
sales <- 200 + (1:36) * 3 + 20 * sin(2 * pi * (1:36) / 12) + rnorm(36, mean = 0, sd = 10)
data <- data.frame(Date = months, Sales = sales)
head(data)## Date Sales
## 1 2020-01-01 207.3952
## 2 2020-02-01 221.0187
## 3 2020-03-01 244.5871
## 4 2020-04-01 230.0256
## 5 2020-05-01 226.2929
## 6 2020-06-01 235.1506Automatic ploting of the hypothetical sales data
ts_data <- ts(data$Sales, start = c(2020, 1), frequency = 12)
autoplot(ts_data, series = "Sales") +
ggtitle("Synthetic Monthly Sales Data") +
xlab("Time") +
ylab("Sales") +
scale_color_manual(values = "blue") + # Customize line color
theme_minimal(base_size = 15) + # Set base font size for better visibility
theme(legend.position = "bottom")
autoplot(object = ts_data,series = "Sales" )
#Check for stationarity
adf_test <- adf.test(ts_data)## Warning in adf.test(ts_data): p-value smaller than printed p-valueprint(adf_test)##
## Augmented Dickey-Fuller Test
##
## data: ts_data
## Dickey-Fuller = -5.3005, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary#Automatically fits a seasonal ARIMA as it sees best.
auto_model <- auto.arima(ts_data)
summary(auto_model)## Series: ts_data
## ARIMA(0,0,0)(1,1,0)[12] with drift
##
## Coefficients:
## sar1 drift
## -0.8392 2.9958
## s.e. 0.0854 0.1095
##
## sigma^2 = 83.1: log likelihood = -93.36
## AIC=192.71 AICc=193.91 BIC=196.25
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4441953 7.126158 4.467821 -0.2139526 1.667343 0.1240441
## ACF1
## Training set 0.1815494We fit SARIMA(1,1,1)(1,1,1)[12] manually.
sarima_model <- Arima(ts_data, order=c(1,1,1), seasonal=c(1,1,1))
summary(sarima_model)## Series: ts_data
## ARIMA(1,1,1)(1,1,1)[12]
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ar1 ma1 sar1 sma1
## 0.0267 -0.7219 -0.8417 -0.0275
## s.e. 0.3068 0.2199 NaN NaN
##
## sigma^2 = 97.05: log likelihood = -91.09
## AIC=192.18 AICc=195.71 BIC=197.86
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3929868 7.156714 4.692756 0.07424912 1.739757 0.1302892
## ACF1
## Training set -0.02246882summary(Arima(ts_data, order=c(1,1,1), seasonal=list(order =c(1,1,1),period =12)))## Series: ts_data
## ARIMA(1,1,1)(1,1,1)[12]
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ar1 ma1 sar1 sma1
## 0.0267 -0.7219 -0.8417 -0.0275
## s.e. 0.3068 0.2199 NaN NaN
##
## sigma^2 = 97.05: log likelihood = -91.09
## AIC=192.18 AICc=195.71 BIC=197.86
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.3929868 7.156714 4.692756 0.07424912 1.739757 0.1302892
## ACF1
## Training set -0.02246882Forecast future values of SARIMA(1,1,1)(1,1,1)[12].
forecasted_values <- forecast(sarima_model, h=12)
autoplot(forecasted_values) +
ggtitle("Sales Forecast for Next 12 Months") +
xlab("Time") +
ylab("Sales") +
theme_minimal()
autoplot(forecasted_values)
Econometrics 03
Simultaneous equation modelling (SEM). Need truffles.csv data for this.
truffles <- read.csv("D:/4th year/AST 431 Statistical Computing X/AST-431 Batch 27/Data/truffles.txt", sep="")
p1_model <- lm(p ~ pf + ps + di,data = truffles)
summary(p1_model)##
## Call:
## lm(formula = p ~ pf + ps + di, data = truffles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.4825 -3.5927 0.2801 4.5326 12.9210
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -32.5124 7.9842 -4.072 0.000387 ***
## pf 1.3539 0.2985 4.536 0.000115 ***
## ps 1.7081 0.3509 4.868 4.76e-05 ***
## di 7.6025 1.7243 4.409 0.000160 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.597 on 26 degrees of freedom
## Multiple R-squared: 0.8887, Adjusted R-squared: 0.8758
## F-statistic: 69.19 on 3 and 26 DF, p-value: 1.597e-12p_hat<- p1_model$fitted.values
p2_model <- lm(q ~ p_hat + ps + di,data = truffles)
summary(p2_model)##
## Call:
## lm(formula = q ~ p_hat + ps + di, data = truffles)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1814 -1.1390 0.2765 1.4595 4.4318
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.27947 3.01383 -1.420 0.167505
## p_hat -0.37446 0.08956 -4.181 0.000291 ***
## ps 1.29603 0.19309 6.712 4.03e-07 ***
## di 5.01398 1.24141 4.039 0.000422 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.68 on 26 degrees of freedom
## Multiple R-squared: 0.6974, Adjusted R-squared: 0.6625
## F-statistic: 19.97 on 3 and 26 DF, p-value: 6.332e-07Dynamic models using built-in data
library(dynlm)## Warning: package 'dynlm' was built under R version 4.4.3library(AER)
data("USMacroG")
plot(USMacroG[, c("dpi", "consumption")],
lty = c(3, 1),
plot.type = "single",
ylab = "")
legend("topleft",
legend = c("income", "consumption"),
lty = c(3, 1), bty = "n")
plot(USMacroG[,c("dpi", "consumption")],plot.type="single",lty=c(2,1))
legend("topleft",legend = c("inco","con"),lty = c(2,1))
Distributed lag models 1st one: consumptiont=β0+β1⋅dpit+β2⋅dpit−1+ut 2nd one: consumptiont=β0+β1⋅dpit+β2⋅dpit−1+β3⋅dpit−2+β4⋅dpit−3+ut
cons_lm1 <- dynlm(consumption ~ dpi + L(dpi), data = USMacroG)
summary(cons_lm1)##
## Time series regression with "ts" data:
## Start = 1950(2), End = 2000(4)
##
## Call:
## dynlm(formula = consumption ~ dpi + L(dpi), data = USMacroG)
##
## Residuals:
## Min 1Q Median 3Q Max
## -190.02 -56.68 1.58 49.91 323.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -81.07959 14.50814 -5.589 7.43e-08 ***
## dpi 0.89117 0.20625 4.321 2.45e-05 ***
## L(dpi) 0.03091 0.20754 0.149 0.882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.58 on 200 degrees of freedom
## Multiple R-squared: 0.9964, Adjusted R-squared: 0.9964
## F-statistic: 2.785e+04 on 2 and 200 DF, p-value: < 2.2e-16cons_lm13 <- dynlm(consumption ~ dpi + L(dpi, 1:3), data = USMacroG)
summary(cons_lm13)##
## Time series regression with "ts" data:
## Start = 1950(4), End = 2000(4)
##
## Call:
## dynlm(formula = consumption ~ dpi + L(dpi, 1:3), data = USMacroG)
##
## Residuals:
## Min 1Q Median 3Q Max
## -199.37 -58.09 1.60 52.32 323.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -85.73584 14.86920 -5.766 3.12e-08 ***
## dpi 0.89301 0.20925 4.268 3.07e-05 ***
## L(dpi, 1:3)1 0.13563 0.28941 0.469 0.640
## L(dpi, 1:3)2 -0.03442 0.28937 -0.119 0.905
## L(dpi, 1:3)3 -0.07246 0.21053 -0.344 0.731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.88 on 196 degrees of freedom
## Multiple R-squared: 0.9964, Adjusted R-squared: 0.9963
## F-statistic: 1.359e+04 on 4 and 196 DF, p-value: < 2.2e-16Autoregressive model: consumptiont=β0+β1⋅dpit+β2⋅consumptiont−1+ut
cons_lm2 <- dynlm(consumption ~ dpi + L(consumption), data = USMacroG)
summary(cons_lm2)##
## Time series regression with "ts" data:
## Start = 1950(2), End = 2000(4)
##
## Call:
## dynlm(formula = consumption ~ dpi + L(consumption), data = USMacroG)
##
## Residuals:
## Min 1Q Median 3Q Max
## -101.303 -9.674 1.141 12.691 45.322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.535216 3.845170 0.139 0.889
## dpi -0.004064 0.016626 -0.244 0.807
## L(consumption) 1.013111 0.018161 55.785 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.52 on 200 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 4.627e+05 on 2 and 200 DF, p-value: < 2.2e-16Granger’s Causality test
# restricted sum of square
g1model <- dynlm(consumption ~ L(consumption), data =USMacroG)
anova(g1model)## Analysis of Variance Table
##
## Response: consumption
## Df Sum Sq Mean Sq F value Pr(>F)
## L(consumption) 1 428664055 428664055 929748 < 2.2e-16 ***
## Residuals 201 92672 461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Unrestricted sum of sqaure
g1model_un <- dynlm(consumption ~ L(consumption) + L(dpi), data =USMacroG)
anova(g1model_un)## Analysis of Variance Table
##
## Response: consumption
## Df Sum Sq Mean Sq F value Pr(>F)
## L(consumption) 1 428664055 428664055 9.5014e+05 < 2e-16 ***
## L(dpi) 1 2440 2440 5.4075e+00 0.02105 *
## Residuals 200 90232 451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Granger causality test
lmtest::grangertest(consumption ~ dpi, data = USMacroG)## Granger causality test
##
## Model 1: consumption ~ Lags(consumption, 1:1) + Lags(dpi, 1:1)
## Model 2: consumption ~ Lags(consumption, 1:1)
## Res.Df Df F Pr(>F)
## 1 200
## 2 201 -1 5.4075 0.02105 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::grangertest(dpi ~ consumption, data = USMacroG)## Granger causality test
##
## Model 1: dpi ~ Lags(dpi, 1:1) + Lags(consumption, 1:1)
## Model 2: dpi ~ Lags(dpi, 1:1)
## Res.Df Df F Pr(>F)
## 1 200
## 2 201 -1 8.327 0.004335 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Econometric VAR model lab e korae nai, but code disilo.
## Vector autoregressive model (VAR)
Var_1 <- dynlm(consumption ~ L(consumption, 1:3) + L(dpi, 1:3), data = USMacroG)
summary(Var_1)##
## Time series regression with "ts" data:
## Start = 1950(4), End = 2000(4)
##
## Call:
## dynlm(formula = consumption ~ L(consumption, 1:3) + L(dpi, 1:3),
## data = USMacroG)
##
## Residuals:
## Min 1Q Median 3Q Max
## -92.724 -10.287 0.822 12.092 44.413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.69704 3.71059 0.188 0.8512
## L(consumption, 1:3)1 1.12289 0.07753 14.484 <2e-16 ***
## L(consumption, 1:3)2 0.10399 0.11124 0.935 0.3510
## L(consumption, 1:3)3 -0.20018 0.08322 -2.405 0.0171 *
## L(dpi, 1:3)1 0.03491 0.05693 0.613 0.5404
## L(dpi, 1:3)2 -0.04252 0.06992 -0.608 0.5438
## L(dpi, 1:3)3 -0.01162 0.05711 -0.203 0.8390
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.51 on 194 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 1.67e+05 on 6 and 194 DF, p-value: < 2.2e-16Var_2 <- dynlm(dpi ~ L(consumption, 1:3) + L(dpi, 1:3), data = USMacroG)
summary(Var_2)##
## Time series regression with "ts" data:
## Start = 1950(4), End = 2000(4)
##
## Call:
## dynlm(formula = dpi ~ L(consumption, 1:3) + L(dpi, 1:3), data = USMacroG)
##
## Residuals:
## Min 1Q Median 3Q Max
## -121.462 -14.712 -2.078 13.750 115.550
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.85925 5.19356 2.091 0.037839 *
## L(consumption, 1:3)1 0.43039 0.10851 3.966 0.000103 ***
## L(consumption, 1:3)2 -0.37263 0.15570 -2.393 0.017653 *
## L(consumption, 1:3)3 -0.02034 0.11648 -0.175 0.861534
## L(dpi, 1:3)1 0.79448 0.07968 9.971 < 2e-16 ***
## L(dpi, 1:3)2 0.24302 0.09786 2.483 0.013862 *
## L(dpi, 1:3)3 -0.06858 0.07994 -0.858 0.391976
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.71 on 194 degrees of freedom
## Multiple R-squared: 0.9997, Adjusted R-squared: 0.9997
## F-statistic: 9.973e+04 on 6 and 194 DF, p-value: < 2.2e-16coeftest(Var_1, vcov. = sandwich)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.697044 2.890480 0.2412 0.80969
## L(consumption, 1:3)1 1.122891 0.091202 12.3122 < 2e-16 ***
## L(consumption, 1:3)2 0.103990 0.111881 0.9295 0.35380
## L(consumption, 1:3)3 -0.200184 0.089856 -2.2278 0.02704 *
## L(dpi, 1:3)1 0.034914 0.056934 0.6132 0.54043
## L(dpi, 1:3)2 -0.042518 0.066909 -0.6355 0.52587
## L(dpi, 1:3)3 -0.011617 0.062566 -0.1857 0.85290
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#install.packages("vars")
library(vars)## Warning: package 'vars' was built under R version 4.4.3## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 4.4.3##
## Attaching package: 'strucchange'## The following object is masked from 'package:stringr':
##
## boundary## Loading required package: urca## Warning: package 'urca' was built under R version 4.4.3library(AER)
data("USMacroG")
var_est <- VAR(as.matrix(cbind(USMacroG[,"consumption"], USMacroG[,"dpi"])), p= 3)
summary(var_est)##
## VAR Estimation Results:
## =========================
## Endogenous variables: USMacroG....consumption.., USMacroG....dpi..
## Deterministic variables: const
## Sample size: 201
## Log Likelihood: -1823.062
## Roots of the characteristic polynomial:
## 1.013 0.9893 0.5221 0.4362 0.3434 0.1721
## Call:
## VAR(y = as.matrix(cbind(USMacroG[, "consumption"], USMacroG[,
## "dpi"])), p = 3)
##
##
## Estimation results for equation USMacroG....consumption..:
## ==========================================================
## USMacroG....consumption.. = USMacroG....consumption...l1 + USMacroG....dpi...l1 + USMacroG....consumption...l2 + USMacroG....dpi...l2 + USMacroG....consumption...l3 + USMacroG....dpi...l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## USMacroG....consumption...l1 1.12289 0.07753 14.484 <2e-16 ***
## USMacroG....dpi...l1 0.03491 0.05693 0.613 0.5404
## USMacroG....consumption...l2 0.10399 0.11124 0.935 0.3510
## USMacroG....dpi...l2 -0.04252 0.06992 -0.608 0.5438
## USMacroG....consumption...l3 -0.20018 0.08322 -2.405 0.0171 *
## USMacroG....dpi...l3 -0.01162 0.05711 -0.203 0.8390
## const 0.69704 3.71059 0.188 0.8512
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 20.51 on 194 degrees of freedom
## Multiple R-Squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 1.67e+05 on 6 and 194 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation USMacroG....dpi..:
## ==================================================
## USMacroG....dpi.. = USMacroG....consumption...l1 + USMacroG....dpi...l1 + USMacroG....consumption...l2 + USMacroG....dpi...l2 + USMacroG....consumption...l3 + USMacroG....dpi...l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## USMacroG....consumption...l1 0.43039 0.10851 3.966 0.000103 ***
## USMacroG....dpi...l1 0.79448 0.07968 9.971 < 2e-16 ***
## USMacroG....consumption...l2 -0.37263 0.15570 -2.393 0.017653 *
## USMacroG....dpi...l2 0.24302 0.09786 2.483 0.013862 *
## USMacroG....consumption...l3 -0.02034 0.11648 -0.175 0.861534
## USMacroG....dpi...l3 -0.06858 0.07994 -0.858 0.391976
## const 10.85925 5.19356 2.091 0.037839 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 28.71 on 194 degrees of freedom
## Multiple R-Squared: 0.9997, Adjusted R-squared: 0.9997
## F-statistic: 9.973e+04 on 6 and 194 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## USMacroG....consumption.. USMacroG....dpi..
## USMacroG....consumption.. 420.6 262.1
## USMacroG....dpi.. 262.1 824.0
##
## Correlation matrix of residuals:
## USMacroG....consumption.. USMacroG....dpi..
## USMacroG....consumption.. 1.0000 0.4451
## USMacroG....dpi.. 0.4451 1.0000(forecast <- predict(var_est))## $USMacroG....consumption..
## fcst lower upper CI
## [1,] 6405.643 6365.446 6445.839 40.19672
## [2,] 6468.504 6407.383 6529.626 61.12183
## [3,] 6535.172 6452.212 6618.131 82.95946
## [4,] 6602.435 6500.487 6704.383 101.94817
## [5,] 6671.212 6551.583 6790.842 119.62967
## [6,] 6740.885 6604.990 6876.781 135.89515
## [7,] 6811.610 6660.420 6962.800 151.18976
## [8,] 6883.271 6717.595 7048.947 165.67580
## [9,] 6955.902 6776.349 7135.454 179.55241
## [10,] 7029.487 6836.543 7222.432 192.94461
##
## $USMacroG....dpi..
## fcst lower upper CI
## [1,] 6688.442 6632.180 6744.704 56.26160
## [2,] 6752.849 6674.420 6831.278 78.42899
## [3,] 6814.348 6715.375 6913.320 98.97210
## [4,] 6879.143 6762.435 6995.852 116.70877
## [5,] 6943.980 6811.522 7076.438 132.45778
## [6,] 7010.201 6863.404 7156.997 146.79648
## [7,] 7077.115 6917.173 7237.057 159.94192
## [8,] 7145.002 6972.805 7317.199 172.19694
## [9,] 7213.728 7030.014 7397.441 183.71392
## [10,] 7283.356 7088.706 7478.005 194.64928# Impulse response
imp <- irf(var_est)
plot(imp)
