# Number of samples
    N <- 10000

    # Draw disturbance variables from normal distributions
    eu <- rnorm(N)*sqrt(3)
    ex <- rnorm(N)
    ey <- rnorm(N)*sqrt(2)
    ez <- rnorm(N)

    # Check them
    hist(eu)
    hist(ex)
    plot(eu,ex)
    
    # Set values of variables
    u <- eu
    x <- -u + ex
    y <- 3*x + ey
    z <- 6*u + y + ez
    
    # Look at pairwise plots
    plot(x,y); 
    plot(y,z);     
    plot(x,z);

    # Fill in matrix B
    B <- matrix(c(0,-1,0,6, 0,0,3,0, 0,0,0,1, 0,0,0,0),4,4)
    print(B)
    
    # Calculate 'total effects' matrix A
    A <- solve(diag(4)-B)
    print(A)

    # Make data frame
    data <- data.frame(u=u, x=x, y=y, z=z)
    
    # Simple regressions, only first one gives the desired causal effect
    model <- lm('y ~ x',data); print(model)
    model <- lm('z ~ y',data); print(model)
    model <- lm('z ~ x',data); print(model)

    # Multiple linear regression, give the desired effects, 
    # but note that u is assumed not observed so these are not really possible!
    model <- lm('z ~ y + u',data); print(model)
    model <- lm('z ~ x + u',data); print(model)

    # Getting the causal effect of y onto z by controlling for x
    model <- lm('z ~ x + y',data); print(model)

    # Covariance matrix of the data
    Cest <- cov(data)
    print(Cest)
    C <- A %*% diag(c(3,1,2,1)) %*% t(A)
    print(C)