v <- cbind(rbind(11,21,31,41),rbind(12,22,32,42))


vec <- function(z) { 
   dimz <- dim(z)
   zret <- t(z[1,])
   for ( i in 2:dimz[1] ) {
     zret <- cbind(zret,t(z[i,])) 
     }
   t(zret)
}

vv <- vec(v)
ii <- diag(x = 1,2,2)

y <- rbind(cbind(0.00454,0.02501),cbind(0.01319,0.01756),cbind(0.02851,-0.00227),
           cbind(0.04018,0.01823))

x <- rbind(cbind(1.0, 0.01673, 0.06252,0.01637),
           cbind(1.0, 0.00794, 0.08857,0.02009),
           cbind(1.0,-0.00188, 0.05403,0.01746),
           cbind(1.0,-0.00877,-0.03892,0.02971))


yy <- vec(y)
xx <- x %x% ii

b <- solve(t(xx) %*% xx) %*% (t(xx)%*% yy)
print(b)

b1 <- lm( y[,1] ~ x[,2:4])
print(b1)
b2 <- lm( y[,2] ~ x[,2:4])
print(b2)

# just for fun

y2 <- vec(t(y))
x2 <- ii %x% x

bb <- solve(t(x2) %*% x2) %*% (t(x2) %*% y2)
print(bb)