library(readxl)
 library(readxl)
 var1 <- read_excel("C:/Users/Tolga OMAY/Desktop/Econometrics Courses/TEKES Dersler/R ile Ekonometri/VAR/var1.xlsx")
 View(var1)
 mydata <- var1
 View(mydata)
### Estimation of the VAR(1) model
 mydata <- var1
 View(mydata)
 
 attach(mydata)
 str(mydata)
 # Defining variables
 y <-ts(mydata[,2],start=c(2000,1),frequency = 12)
 x <-ts(mydata[,3],start=c(2000,1),frequency = 12)
 plot(y)
 TT <- dim(mydata)
 TT
Var_1 <- ar(cbind(x, y), method="ols", dmean=T, intercept=F)
Var_1$ar

data <- cbind(x, y)
View(data)
#Ordinary least squares is used to fit the model to the mean adjusted series – with
#dmean set to TRUE and intercept set to FALSE since the latter parameter will
#not be required.

acf(Var_1$res[-c(1:1), 1])  
acf(Var_1$res[-c(1:1), 2])

# lag 1 cıktı 1:1 3 cıksa :3
#From the code above, we see that the best-fitting VAR model is of order 1.
#The correlogram of the residual series indicates that the residuals are approximately
#install.packages("vars")
library(vars)
#?VAR
VAR1.var <- VAR(cbind(x, y),p=1, type = "both", ic = c("HQ"))
coef(VAR1.var)
#bivariate white noise, thus validating the assumptions for a VAR model

acf(resid(VAR1.var)[, 1])
acf(resid(VAR1.var)[, 2])


Var1.pred <- predict(VAR1.var, n.ahead = 4)
Var1.pred

y.pred <- ts(Var1.pred$fcst$y[, 1], st = 1988, fr = 4)
x.pred <- ts(Var1.pred$fcst$x[, 1], st = 1988, fr = 4)
ts.plot(cbind(window(y, start = 1981), y.pred), lty = 1:2)
ts.plot(cbind(window(x, start = 1981), x.pred), lty = 1:2)

gdp <-ts(y, start=c(2008:1),frequency=12)
unemp <-ts(x, start=c(2008:1),frequency=12)

#####using more library VAR


library(urca)
#install.packages("mFilter")
library(mFilter)
#install.packages("TSstudio")
library(TSstudio)
library(tidyverse)
library(tseries)
library(forecast)
plot(y, x)
autoplot(cbind(gdp,unemp))
VAR1.mod <- cbind(gdp,unemp)
colnames(VAR1.mod) <- cbind("DGP","Unemployment") 
lagselect <- VARselect(VAR1.mod, lag.max= 10, type="const")
lagselect$selection
VAR1.mod1 <- VAR(VAR1.mod, p = 1, type = "const", season = NULL, exog =NULL)
summary(VAR1.mod1)


#Diagnosis of the VAR

# serial correlation
serial1 <- serial.test(VAR1.mod1, lags.pt = 12 , type = "PT.asymptotic")
serial1
# there is no serial correlation p value is 0.3117

# heteroscedasticity
 ARCH1 <- arch.test(VAR1.mod1, lags.multi = 12 , multivariate.only = TRUE)
ARCH1
# there is no heteroscedasticity p value is 0.3211

#Normality Test

Norm1 <- normality.test(VAR1.mod1,multivariate.only = TRUE )
Norm1

# testing for Structural or model stability test

Stability1 <- stability(VAR1.mod1, type = "OLS-CUSUM")
plot(Stability1)


# Granger Causality

GrangerGDP <- causality(VAR1.mod1, cause = "DGP")
GrangerGDP

GrangerUnemp <- causality(VAR1.mod1, cause = "Unemployment")
GrangerUnemp

# Impulse Response Function

GDPirf  <- irf(VAR1.mod1,impulse ="Unemployment", response= "DGP", n.ahead = 20, boot = TRUE)
plot(GDPirf, ylab ="DGP", main = "Shock Unemployment")

Unempirf  <- irf(VAR1.mod1,impulse ="DGP", response= "Unemployment", n.ahead = 20, boot = TRUE)
plot(Unempirf, ylab ="Unemployment", main = "Shock GDP")


# Variance Decompoistion

FEVDI <- fevd(VAR1.mod1, n.ahead= 10)
plot(FEVDI)

#¦ VAR forecasting

forecast <- predict(VAR1.mod1, n.ahead= 20, ci = 0.95)
fanchart(forecast, names = "DGP")
fanchart(forecast, names = "Unemployment", colors = "Pink")
#?fanchart


#1 Structural VAR
set.seed(0)
y<-arima.sim(model=list(ar=.75), n=120, innov = runif(120, 0, 1))
z<-y-0.5

par(mfrow=c(1,1))
ts.plot(z)
inf <-ts(z, start=c(2008:1),frequency=12)

onat <- diag(3)
onat
onat[2,1] <- NA
onat[3,1] <- NA
onat[3,2] <- NA
onat
# sırayla y ikisini etkiliyor, x  

sv <- cbind(gdp,unemp,inf)
colnames(sv) <- cbind("DGP", "Unemployment", "inflation")
View(sv)
autoplot(sv)
# lagselect
lagselect1 <- VARselect(sv, lag.max = 4, type = "const")
lagselect1 <- VARselect(sv, lag.max = 4, type = "both")
lagselect1 <- VARselect(sv, lag.max = 4, type = "none", season = NULL)
lagselect1$selection
Model1 <- VAR(sv, p=1, season = NULL, exog = NULL, type = "const")
SVARMod1 <- SVAR(Model1, Amat = onat, Bmat = NULL, hessian = TRUE, estmethod = c("scoring","direct"))
SVARMod1

#SVAR implres

SVARirf  <- irf(SVARMod1,impulse ="Unemployment", response= "DGP")
plot(SVARirf, ylab ="DGP", main = "Shock Unemployment")
SVARirf  <- irf(SVARMod1,impulse ="Unemployment", response= "inflation")
plot(SVARirf, ylab ="inflation", main = "Shock Unemployment")
SVARirf  <- irf(SVARMod1,impulse ="DGP", response= "Unemployment")
plot(SVARirf, ylab ="Unemployment", main = "Shock DGP")
SVARirf  <- irf(SVARMod1,impulse ="DGP", response= "inflation")
plot(SVARirf, ylab ="inflation", main = "Shock DGP")
SVARirf  <- irf(SVARMod1,impulse ="DGP", response= "DGP")
plot(SVARirf, ylab ="DGP", main = "Shock DGP")