Granger on vector error-correction model Say I have an vector error-corretion model for $3$ time series $x,y,z$.
$$x_t=\alpha_1 +\sum_{i=1}^{p}\theta_{1i} y_{t-i}+\sum_{i=1}^{p}\delta_{1i} z_{t-i}+\beta_1 ECT_{t-1}+e_{1t}$$
$$y_t=\alpha_2 +\sum_{i=1}^{p}\theta_{2i} y_{t-i}+\sum_{i=1}^{p}\delta_{2i} z_{t-i}+\beta_2 ECT_{t-1}+e_{2t}$$
$$z_t=\alpha_3 +\sum_{i=1}^{p}\theta_{3i} y_{t-i}+\sum_{i=1}^{p}\delta_{3i} z_{t-i}+\beta_3 ECT_{t-1}+e_{3t}$$
where $ETC$ denotes error term and $e_{it}$ is the stochastic error term with mean zero and constant variance.
I am going to run Granger causality test on it in R. But I not fully understand should I estimate $\alpha, \theta, \delta $ coefficients before Granger causality test or Granger causality test can receive variable coefficients. 
How can I apply the Granger causality test on these equations, this case?
 A: You should test for Granger (non-)casuality in the underlying VAR-model in (log) levels, rather than the VECM representation of it. This code will reproduce (part of) the example from Dave Giles' blog post using lm() and the lmtest package. I found the data here (where another implementation in R is also suggested)
library(lmtest)
library(xts)    

### READ AND FORMAT DATA ###

coffee <- read.csv("http://christophpfeiffer.org/wp-content/uploads/2012/11/coffee_data.csv", header=T,sep=";")
coffee <- coffee[1:615,]
colnames(coffee) <- c("time", "arabica", "robusta")
coffee$time <- as.yearmon(coffee$time, "%YM%m")
coffee <- coffee[coffee$time>=as.yearmon("1975M06", "%YM%m"),]
coffee <- xts(coffee[,-1], order.by = coffee$time)

### MAKE LAGS ###

coffee.l7 <- lag(coffee, 7)
colnames(coffee.l7) <- paste(colnames(coffee.l7), ".l7", sep = "")
arabica.l1l6 <- lag(coffee[,"arabica"], 1:6)
colnames(arabica.l1l6) <- paste("arabica.l", 1:6, sep = "")
robusta.l1l6 <- lag(coffee[,"robusta"], 1:6)
colnames(robusta.l1l6) <- paste("robusta.l", 1:6, sep = "")

### TODA-YAMAMOTO GRANGER-TEST, H0: ARABICA IS A NON-CAUSE OF ROBUSTA ###

unrestricted.fit <- lm(coffee[,"robusta"] ~ arabica.l1l6 + robusta.l1l6 + coffee.l7)
restricted.fit <- lm(coffee[,"robusta"] ~ robusta.l1l6 + coffee.l7)
waldtest(restricted.fit, unrestricted.fit, test = "Chisq")

Extending it to your case, with more than two variables, I would guess that we would do it like this; with a VAR(2) for labour productivity, real wage, employment, and unemployment rate and assuming all variables ar I(1):
library(vars)
data(Canada)
Canada <- as.xts(Canada)

### MAKE LAGS ###

prod.l1l2 <- lag(Canada[,"prod"], 1:2)
colnames(prod.l1l2) <- paste("prod.l", 1:2, sep = "")
rw.l1l2 <- lag(Canada[,"rw"], 1:2)
colnames(rw.l1l2) <- paste("rw.l", 1:2, sep = "")
U.l1l2 <- lag(Canada[,"U"], 1:2)
colnames(U.l1l2) <- paste("U.l", 1:2, sep = "")
e.l1l2 <- lag(Canada[,"e"], 1:2)
colnames(e.l1l2) <- paste("e.l", 1:2, sep = "")
Canada.l3 <- lag(Canada, 3)
colnames(Canada.l3) <- paste(colnames(Canada.l3), "l3", sep = ".")

### TODA-YAMAMOTO GRANGER-TEST, H0: LABOUR PRODUCTIVITY IS A NON-CAUSE OF REAL WAGE ###

ur.fit <- lm(Canada[,"rw"] ~ rw.l1l2 + prod.l1l2 + U.l1l2 + e.l1l2 + Canada.l3)
r.fit <- lm(Canada[,"rw"] ~ rw.l1l2 + U.l1l2 + e.l1l2 + Canada.l3)
waldtest(r.fit, ur.fit, test = "Chisq")

Since it is not very practical to lag the data manually, you could use VAR() in the vars package instead of lm(), which also includes the ca.jo() function from urca for VECM and Johansen's test. You can confirm that VAR() fits the same model as lm() above. causality(), however, only reports the F-version of the test (so it was not able to check it against Giles' results); also, I am not sure how it handles the exogenous variables from the exogen argument.
var.fit <- VAR(coffee, type = "const", p = 6)
var.fit_alt <- VAR(coffee, type = "const", p = 6, exogen = coffee.l7) # INCLUDE LAG 7 AS EXOGENOUS VARIABLES FOR TY-TEST

cbind(vars = var.fit_alt$varresult$robusta$coefficients,
      lm = unrestricted.fit$coefficients[c(2,8,3,9,4,10,5,11,6,12,7,13,1,14,15)])

causality(var.fit_alt, cause = "arabica")$Granger

