library(OpenMx)
mxOption(key='Number of Threads', value=parallel::detectCores() - 1)

# arrow = 1  from -> to
# arrow = 2  from <-> to
# free = FALSE -> fixed path (parameter values are ex ante fix)
# free = TRUE is a free path (parameter values are estimated)

# from: sources of the new paths
# values: Starting values of parameters
# connect: type of source to sink connection
# labels: name of the paths

# values is the starting value of a path # reasonable values for faster convergence

# first try to replicate model in
# "Linear dynamic panel-data estimation using maximum likelihood and structural equation modeling"

# yit = λyit−1 + x∗itβ + w∗iδ + αi + ξt + υit ∀i, t
# E(υit|yt−1,i , xt,i, wi, αi) = 0 ∀i, 


library(haven)
wages <- haven::read_dta("https://www3.nd.edu/~rwilliam/statafiles/wageswide.dta")
#is(wages)
wages <- as.data.frame(scale(wages))
#wages <- as.matrix(wages)
#is(wages)
#View(wages)

head(wages)
summary(wages)

## 595 households who reported a non-zero wage (measured for 7 years: 1976 - 1982)
# y = wks: number of weeks employed in each year 
# x = union (dummy): 1 if wage set by union contract
# w = lwage: ln of wage in each year
# z = ED: years of education in 1976


manifests <- c("wks2", "wks1", "ed", "lwage1", "union1")

latents <- "alpha"

cov_data <- cov(wages[,manifests])


model1 <- mxModel("Test1", type = "RAM",
                  
                  # Manifest variables
                  manifestVars = manifests,
                  
                  # Latent Variables
                  latentVars = latents,
                  
                  mxPath(from = latents,
                         to = "wks2",
                         arrows = 1,
                         free = FALSE, 
                         values = 1),
                  
                  
                  mxPath(from = c("ed", "wks1", "lwage1", "union1"),
                         to = "wks2", 
                         arrows = 1,
                         free = TRUE, values = rep(1, 4)),
                  
                  mxPath(from = "wks2",
                         to = "wks2",
                         arrows = 2, connect = "single",
                         free = TRUE, values = 1), # 26.366 according to mxAutoStart(model1) -> below
                  
                  #mxPath(from = latents,
                  #       to = c("wks1", "lwage1", "union1"),
                  #       arrows = 2, connect = "single",
                  #       free = FALSE, values = c(1,1,1)),
                  
                  #mxPath(from = c("ed", "wks1", "lwage1", "union1"),
                  #       arrows = 2, connect = "unique.pairs",
                  #       free = TRUE, values = c(1,1,1,1)),
                  
                  mxPath(from = "one",
                         to = "wks2",
                         values = 1),
                  
                  mxData(wages, type = "raw", numObs = 595))

model1

round(eigen(mxGetExpected(model1, "covariance"))$values, 3)

smodel1 <- mxAutoStart(model1)
mxGetExpected(smodel1, "covariance")



fitmodel <- mxRun(model1)
summary(fitmodel)