Synthetique v0.Rmd

```{r}
library(ggplot2)
library(reshape2)
library(dplyr)
library("Synth")
library(tidyverse)
library(numDeriv)
library('alabama')
library("rgenoud")
library(tidyr)
library(optimx)
library(kernlab)

```

```{r}
capital=read.csv("Statistical_Capital.csv")
croiss=read.csv("Statistical_Growth-Accounts.csv")
travail=read.csv("Statistical_Labour.csv")

compte=read.csv("Statistical_National-Accounts.csv")
description=read.csv("Variable-Description.csv")

autre_pays=read.csv("WB_data2.csv")
pays2lettres=read.csv('liste_pays_2lettres.csv')
pays3lettres=read.csv('liste_pays_3lettres.csv')
codepays=read.csv("codes_pays.csv")
taux_change <- read_csv("taux_change.csv")
donnee=read.csv('donnee.csv')
```

```{r}
donnees_bis=read.csv('donnees_complementaires.csv')
```

#les lignes avec des données
```{r}
lignedata=c('X1995', 'X1996', 'X1997', 'X1998', 'X1999', 'X2000', 'X2001', 'X2002', 'X2003', 'X2004', 'X2005', 'X2006', 'X2007', 'X2008', 'X2009', 'X2010', 'X2011', 'X2012', 'X2013', 'X2014', 'X2015', 'X2016', 'X2017')#les colonnes de données 
date=seq(2000, 2017)
```


```{r}
codepays=codepays %>% rename( country=ISO2.digit.Alpha)

compte=compte %>% mutate(country = str_replace_all(compte$country, c("EL" = "GR", "UK" = "GB"))) 
comptebis=left_join(compte,codepays,by="country")
compte=comptebis%>% rename(country3L=ISO3.digit.Alpha)

capital=capital %>% mutate(country = str_replace_all(capital$country, c("EL" = "GR", "UK" = "GB"))) 
capitalbis=left_join(capital,codepays,by="country")
capital=capitalbis%>% rename(country3L=ISO3.digit.Alpha)

croiss=croiss%>% mutate(country = str_replace_all(croiss$country, c("EL" = "GR", "UK" = "GB")))
croissbis=left_join(croiss,codepays,by="country")
croiss=croissbis%>% rename(country3L=ISO3.digit.Alpha)

```
```{r}
donnees_bis=donnees_bis%>% rename(taux_imposition=taux.d.imposition.sur.bÃ.nÃ.fices)
donnees_bis=donnees_bis%>% rename(salaire_min=salaire.horaire.minimum.en.US..PPA)
donnees_bis=donnees_bis%>% rename(protection_emploi=indice.protection.de.l.emploi)
donnees_bis=donnees_bis%>% rename(diplome=X..des.25.64ans.diplomes.du.superieur)
donnees_bis=donnees_bis%>% rename(IDE=IDE.entrants.en...PIB)
```

#fonction qui permet de récupérer des dataframes utilisables sur une donnée
```{r}
melt_ada <- function(vari,codei,basei=compte,tps=lignedata[3:22]){
tableau=subset(basei,(var==vari)& code==codei)[c("country3L", tps,"var")]
indligneNA <- which(is.na(tableau),arr.ind=TRUE)[,1]
tableau=tableau[-indligneNA,]
df.melted <- melt(tableau,id=c('country3L','var'))
df.melted$temps=as.numeric(str_sub(df.melted$variable, 2, 5))
df.melted$country=as.character(df.melted$country3L)
return (df.melted[-3]%>% select(-country3L))}
```


```{r}
df.melted=inner_join(melt_ada('VA','TOT'),melt_ada("EMP",'TOT'),by=c("country","temps"))
compte%>%filter(country=="JP",code=="TOT",var=="EMP") #on perd le japon a cause de 2 NA sur l'emploi

df.melted$index=c(1:length((unique(df.melted$country))))
mode(df.melted[,'index'])
df.melted$value.x=as.numeric(df.melted$value.x)#on transtype pour permettre aux fonctions de tourner
df.melted$country=as.character(df.melted$country)

df.melted=inner_join(df.melted, melt_ada("LP_QI","TOT"), by=c("country","temps"))
df.melted=inner_join(df.melted, melt_ada("LAB","TOT",croiss), by=c("country","temps"))

df.melted=df.melted%>% rename(VA=value.x,Emploi=value.y,Productivite=value.x.x,LAB=value.y.y)
```

#le calcul du salaire réel horaire
```{r}
df.W=inner_join(melt_ada('COMP','TOT'), melt_ada("VA_PI",'TOT'),by=c("country","temps"))
df.W=inner_join(df.W, melt_ada("H_EMPE",'TOT'),by=c("country","temps"))
df.W$index=c(1:length((unique(df.W$country))))
df.W$country=as.character(df.W$country)
df.W$W=(df.W$value.x*100/df.W$value.y)*1000/df.W$value#le nombre d'heures travaillées est en milliers
df.W$VA=as.numeric(df.W$W)
df.W=df.W[,c('country','W', 'temps')]

df.melted2=inner_join(df.melted, df.W, by=c("country","temps"))
df.melted2$index=c(1:length((unique(df.melted2$country))))
```


#Intégration de la population taux de chomage et exportations

```{r}
df.melted2=df.melted2[c("country",'temps','index','VA','Emploi','W','Productivite','LAB')]
dfpop=subset(autre_pays,Simple.Indicator=="pop")[c("Country.Code","Simple.Indicator",lignedata)]
rownames(dfpop)<-dfpop[,1]

df_pop<-t(dfpop[-c(1,2)])
df <- data.frame(x=1995:2017, data.frame(df_pop))
df.pop <- melt(df, id="x")%>%rename(country=variable, temps=x,pop=value)
```


```{r}
dfchom=subset(autre_pays,Simple.Indicator=="unemployment rate")[c("Country.Code","Simple.Indicator",lignedata)]
rownames(dfchom)<-dfchom[,1]
dfexports=subset(autre_pays,Simple.Indicator=="exports_hab")[c("Country.Code","Simple.Indicator",lignedata)]
rownames(dfexports)<-dfexports[,1]

df_chom<-t(dfchom[-c(1,2)])
df <- data.frame(x=1995:2017, data.frame(df_chom))
df.chom <- melt(df, id="x")%>%rename(country=variable, temps=x,chomage=value)
df_exp<-t(dfexports[-c(1,2)])
df <- data.frame(x=1995:2017, data.frame(df_exp))
df.exp <- melt(df, id="x")%>%rename(country=variable, temps=x,export=value)

df.melted2=left_join(df.melted2,df.chom)
df.melted2=left_join(df.melted2,df.pop)
df.melted2=left_join(df.melted2,df.exp)

df.melted2$VA_pop=df.melted2$VA/df.melted2$pop*1000000
df.melted2$Part_lab=df.melted2$LAB/df.melted2$VA
```

#Investissement
```{r}
df.invest=melt_ada('Iq_GFCF','TOT',capital,tps=lignedata[1:22])
df.melted_i=inner_join(df.melted2,df.invest,by=c("country","temps"))
df.melted_i$index=c(1:length((unique(df.melted_i$country))))
```

#Donnée réelle et indice 100
```{r}
taux_change2010 <- taux_change %>% filter(year == 2010) %>%
  rename(`rate2010 in NAT/US$` = `rate in NAT/US$`, country=country3) %>% select(-year)

df.melted_i <- df.melted_i  %>% inner_join(taux_change2010, by = "country")
df.melted_i<-df.melted_i %>% mutate(VA_pop = VA_pop/(pop*`rate2010 in NAT/US$`))

donn2010=df.melted_i%>%filter(temps==2010)%>%select('country','VA','Emploi')%>%rename(VA2010=VA, Emp2010=Emploi)
df.melted_i=df.melted_i%>%left_join(donn2010,on="country")

df.melted_i$VA_ind=(df.melted_i$VA/df.melted_i$VA2010)*100
df.melted_i$Emp_ind=(df.melted_i$Emploi/df.melted_i$Emp2010)*100

df.melted_i$Invest_pop=df.melted_i$value/df.melted_i$pop
```


```{r}
donnees_bis
df.fisc=inner_join(df.melted_i,donnees_bis%>%rename(country=country3, temps=year))
indligneNA <- which(is.na(df.fisc),arr.ind=TRUE)[,1]
df.fisc[unique(indligneNA),][c(1:2,22:ncol(df.fisc))]%>%arrange(by=temps)
df.fisc=df.fisc[df.fisc$temps>1999,]#on a pas tous les pays pourles premières années, donc on prend à partir de 2000
df.fisc$index=c(1:length((unique(df.fisc$country))))
```

#Evolution de la VA : base df.evol qui commence en 98 à cause du lag.

```{r}
f=function(variabl=c('VA'), df.melt=df.melted_i, df.evol=df.evol){
base=NULL
#On peut facilement rajouter les évolutions des autres var en changeant un tt petit peu le code

for (pays in unique(df.melt$country)){
a=df.melt%>%filter(country==pays)%>%select(temps,country, variabl)

a[,variabl][2:nrow(a)]=diff(a[,variabl])/a[,variabl][2:nrow(a)]*100
a=a%>%rename(evol_VA=variabl)
base=rbind(base,a)}
df.evol=df.melt%>%full_join(base,on=c('country','temps'))
df.evol=df.evol%>%filter(temps>1997)
return(df.evol)}
df.evol=f('Emploi')
df.evol=df.evol%>%rename(evol_Emp=evol_VA)
df.evol=f('VA', df.melt=df.evol)
```

###Le calcul des poids optimaux pour l'ensemble de la période d'optimisation


Préparation des donnnées pour l'optimisation
```{r}
quiet <- function(x) { 
  sink(tempfile()) 
  on.exit(sink()) 
  invisible(force(x)) 
} 

```

```{r}
calculpoids=function(df.melted, value, pays='FRA', poids_max=TRUE){
  #préparation des données
pays_autres=length(c(table(df.melted$country)))-1
Xt=subset(df.melted, country==pays& temps<2013 &temps>1999)
Xt=Xt[,value]#on récupère les données pour la France

Xn=subset(df.melted, temps<2013 &temps>1999)
Xn=Xn[Xn$country!=pays,]#on récupère les autres
Xn=matrix(Xn[,value], nrow=pays_autres)

#la fonction d'optimisation

x0=c(rep(1/pays_autres,pays_autres))#on crée un vecteur de départ au hasard
f=function(x){norm(Xt-x%*%Xn, type="2")}
n=function(x){x}#contrainte de positivité des poids (hin(x)>0)
if(poids_max==TRUE){n=function(x){c(x, 0.25-x)}}
q=function(x){sum(x)-1}#fonction de la contrainte d'égalité (heq(x)=0)
opti=quiet(auglag(x0, f, hin=n, heq=q))#fonction d'optimisation sur toute la période
opti$par=matrix(opti$par, ncol=1)}
```

```{r}
gen_liste=function(liste_value, date=c(list(2000:2007), list(2008:2011), 2012))
{
a=0

for(i in liste_value ){
  for (j in date){
    a=a+1
    if(a==1){liste=list(list(i,j, "mean"))}
    
    else{liste=c(liste, list(list(i, j, "mean")))  }
  }
}
liste=as.list(liste)

return(liste)
}

```
###la fonction avec stratification du synthétique

```{r}
df.melt=df.melted_i
liste_value='VA_ind'
pays='FRA'
date0=1997
date1=2005
date2=2012
date3=2016
synthdon=TRUE
synth.out=0
  valeur=liste_value[1]
  #on récupère les poids
  ind_pays=as.numeric(df.melt[df.melt$country==pays,]$index[1])
  if(synthdon){dataprep.out=Synth::dataprep(foo=df.melt,
                predictors = liste_value, #Les variables utilisées en contrôle
                predictors.op = "mean", #La manière dont on les prend en compte
                time.predictors.prior = 2012, #période pré-traitement
                dependent=liste_value[1], #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = ind_pays, # n° de l'unité traitée
                controls.identifier = c(1:max(df.melt$index))[-ind_pays], #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = date1:date2, #période de maximisation
                time.plot=date0:date3) #période tracée sur le graphique
synth.out <- Synth::synth(data.prep.obj = dataprep.out, method = "BFGS")}else{synth.out=synth.out}
```

```{r}

solpoids=synth.out$solution.w
#puis on construit les séries augmentées
a=1


for(poids in solpoids){
if(a!=ind_pays){

  newrow=data.frame(valeur=df.melt[df.melt$index==a & df.melt$temps==2016,][,valeur]*poids, index=a,country=unique(df.melt[df.melt$index==a,]$country))
  if(a==1){dfstrat=newrow} 
  else{dfstrat=rbind(dfstrat, newrow)}}
a=a+1
}

numgrand=order(dfstrat$valeur, decreasing=TRUE)
numbis=dfstrat$index[numgrand]

#dfstrat=dfstrat[order(dfstrat$valeur, decreasing=TRUE),]

#for (num in 1:length(numgrand)){
  #if(numgrand[num]>=ind_pays){numgrand[num]=numgrand[num]+1}
  #}#on compense l'absence de la France

newrow=df.melt[df.melt$index==dfstrat$index[numgrand[1]] & df.melt$temps>date0-1 & df.melt$temps<date3+1,][,valeur]*solpoids[dfstrat$index[numgrand[1]]]
newdata=data.frame(valeur=newrow, country=unique(dfstrat[dfstrat$index==dfstrat$index[numgrand[1]],]$country), temps=date0:date3)
dfstratfin=newdata

for(pays in 2:length(numgrand)){#(max(df.melt$index)-2)
   newrow=newrow+df.melt[df.melt$index==dfstrat$index[numgrand[pays]]& df.melt$temps>date0-1 & df.melt$temps<date3+1,][,valeur]*solpoids[dfstrat$index[numgrand[pays]]]
   if (pays<6){
     newdata=data.frame(valeur=newrow, country=unique(dfstrat[dfstrat$index==dfstrat$index[numgrand[pays]],]$country), temps=date0:date3)
   
     dfstratfin=rbind(dfstratfin, newdata)
   }
   if(pays==max(df.melt$index)-2){
     newdata=data.frame(valeur=newrow, country='autres', temps=date0:date3)
     dfstratfin=rbind(dfstratfin, newdata)
     }
}
newrow=df.melt[df.melt$index==ind_pays& df.melt$temps>date0-1 & df.melt$temps<date3+1,][,valeur]
newdata=data.frame(valeur=newrow, country='Vraie valeur', temps=date0:date3)
dfstratfin=rbind(dfstratfin, newdata)
ggplot(data=dfstratfin,aes(x=temps, valeur)) +
  geom_line(aes(group=country, col=country))  

```
```{r}
synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
poids1=synth.tables$tab.w%>%rename(poids=w.weights)
poids1 = poids1%>%rename(country=unit.names)
poids1 = poids1[order(poids1$poids, decreasing=TRUE),]
```

```{r}
df.melt=df.melted_i
liste_value='VA_ind'
pays='FRA'
date0=2000
date1=2005
date2=2012
date3=2016
ordonnee = NA
nb_country_plot = 5
synthdon=TRUE
synth.out=synth.out
dataprep.out=dataprep.out
valeur=liste_value[1]
ind_pays= unique(subset(df.melt, country == pays)['index'])[1, 1]
dateopti = date0:date3

if(synthdon==FALSE){
dataprep.out=Synth::dataprep(foo=df.melt,
                predictors.op = "mean", #La manière dont on les prend en compte
                time.predictors.prior = date2, #période pré-traitement
                special.predictors = gen_liste(liste_value, c(date1:2009,2011:date2)),
                dependent=liste_value[1], #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = ind_pays, # n° de l'unité traitée
                controls.identifier = c(unique(subset(df.melt, country!= pays)['index'])[,1]), #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = dateopti, #période de maximisation
                time.plot=date0:date3) #période tracée sur le graphique

synth.out <- Synth::synth(data.prep.obj = dataprep.out, method = "BFGS")


Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = liste_value[1], Xlab = "year", Legend = c(pays, paste('synthetic',pays)), Legend.position = "bottomright", Ylim = ordonnee)}
synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
poids1=synth.tables$tab.w%>%rename(poids=w.weights)
poids1 = poids1%>%rename(country=unit.names)
poids1 = poids1[order(poids1$poids, decreasing=TRUE),]
solpoids=synth.out$solution.w
dfstrat = data.frame()

for (countries in unique(poids1$country)){
        for (i in 1:length(unique(df.melt$country))){ 
          #On sépare maintenant selon si le pays a un indice plus petit que celui d'intérêt, ou plus grand, pour bien attribuer le poids correspondant au bon indice 
          
          if (countries == df.melt$country[i] & i <ind_pays){
            for (years in date0:date3){
            newrow = data.frame(valeur=df.melt[df.melt$country==countries & df.melt$temps==years,][,valeur]*solpoids[i], country=countries, index=unique(df.melt[df.melt$country==countries,]$index), weight = solpoids[i], temps = years)
            dfstrat = rbind(dfstrat, newrow)}}
          
          if (countries == df.melt$country[i] & i >ind_pays){
            for (years in date0:date3){
            newrow = data.frame(valeur=df.melt[df.melt$country==countries & df.melt$temps==years,][,valeur]*solpoids[i-1],country=countries, index=unique(df.melt[df.melt$country==countries,]$index),  weight = solpoids[i-1], temps = years)
            dfstrat = rbind(dfstrat, newrow)}}
          }  
}


#On crée maintenant le df avec les valeurs cumulées des contributions de chaque pays
dfstratfin = data.frame()
for (years in date0:date3){#on initialise
  newrow = dfstrat[dfstrat$country == unique(dfstrat$country)[1] & dfstrat$temps == years,]$valeur
  temp_df = data.frame(valeur=newrow, country=unique(dfstrat$country)[1], temps = years, index=unique(df.melt[df.melt$country==unique(dfstrat$country)[1],]$index), 
weight = dfstrat[dfstrat$country == unique(dfstrat$country)[1] & dfstrat$temps == years, ]$weight)
  dfstratfin = rbind(dfstratfin, temp_df)
  
  for (countries in unique(dfstrat$country)[2:nb_country_plot]){ #on garde ceux les plus élevés
    newrow = newrow + dfstrat[dfstrat$country == countries & dfstrat$temps == years, ]$valeur
    
    temp_df = data.frame(valeur=newrow, index=unique(df.melt[df.melt$country==countries,]$index), country=countries, 
weight = dfstrat[dfstrat$country == countries & dfstrat$temps == years,]$weight, temps = years)
    dfstratfin = rbind(dfstratfin, temp_df)
  }}



#On rajoute l'unité synthétique
for (years in date0:date3){
  newrow = dfstrat[dfstrat$country == unique(dfstrat$country)[1] & dfstrat$temps == years,]$valeur
  for (countries in unique(dfstrat$country)[2:length(unique(dfstrat$country))]){ #on rajoute les petites valeurs
    newrow = newrow + dfstrat[dfstrat$country == countries & dfstrat$temps == years, ]$valeur}
  temp_df = data.frame(valeur=newrow, index=0, country='Synthetic Unit', weight = 1, temps = years)
  dfstratfin = rbind(dfstratfin, temp_df)
}
dfstratfin = dfstratfin[order(dfstratfin$weight, decreasing=TRUE),]


#On ajoute la vraie valeur pour la France
for (years in date0:date3){
  temp_df = data.frame(valeur=df.melt[df.melt$country == pays & df.melt$temps == years,][,valeur], country= pays, temps=years, weight = 1, index = ind_pays)

  dfstratfin=rbind(dfstratfin, temp_df)}

colorplot = ggplot(data=dfstratfin,aes(x=temps, y=valeur))+
  geom_line(aes(group=country, col=country))

print(colorplot)
```
```{r}
df_stratTitouan()
```

```{r}
df_stratTitouan=function(df.melt=df.melted_i, liste_value='VA_ind', pays='FRA', date0=2000, date1=2005, date2=2012, date3=2016, ordonnee = NA, nb_country_plot = 5, synthdon=FALSE, synth.out=0, dataprep.out=0){
valeur=liste_value[1]
ind_pays= unique(subset(df.melt, country == pays)['index'])[1, 1]
dateopti = date0:date3

if(synthdon==FALSE){
dataprep.out=Synth::dataprep(foo=df.melt,
                predictors.op = "mean", #La manière dont on les prend en compte
                time.predictors.prior = date2, #période pré-traitement
                special.predictors = gen_liste(liste_value, c(date1:2009,2011:date2)),
                dependent=liste_value[1], #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = ind_pays, # n° de l'unité traitée
                controls.identifier = c(unique(subset(df.melt, country!= pays)['index'])[,1]), #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = dateopti, #période de maximisation
                time.plot=date0:date3) #période tracée sur le graphique

synth.out <- Synth::synth(data.prep.obj = dataprep.out, method = "BFGS")


Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = liste_value[1], Xlab = "year", Legend = c(pays, paste('synthetic',pays)), Legend.position = "bottomright", Ylim = ordonnee)}
synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
poids1=synth.tables$tab.w%>%rename(poids=w.weights)
poids1 = poids1%>%rename(country=unit.names)
poids1 = poids1[order(poids1$poids, decreasing=TRUE),]
solpoids=synth.out$solution.w
dfstrat = data.frame()

for (countries in unique(poids1$country)){
        for (i in 1:length(unique(df.melt$country))){ 
          #On sépare maintenant selon si le pays a un indice plus petit que celui d'intérêt, ou plus grand, pour bien attribuer le poids correspondant au bon indice 
          
          if (countries == df.melt$country[i] & i <ind_pays){
            for (years in date0:date3){
            newrow = data.frame(valeur=df.melt[df.melt$country==countries & df.melt$temps==years,][,valeur]*solpoids[i], country=countries, index=unique(df.melt[df.melt$country==countries,]$index), weight = solpoids[i], temps = years)
            dfstrat = rbind(dfstrat, newrow)}}
          
          if (countries == df.melt$country[i] & i >ind_pays){
            for (years in date0:date3){
            newrow = data.frame(valeur=df.melt[df.melt$country==countries & df.melt$temps==years,][,valeur]*solpoids[i-1],country=countries, index=unique(df.melt[df.melt$country==countries,]$index),  weight = solpoids[i-1], temps = years)
            dfstrat = rbind(dfstrat, newrow)}}
          }  
}


#On crée maintenant le df avec les valeurs cumulées des contributions de chaque pays
dfstratfin = data.frame()
for (years in date0:date3){#on initialise
  newrow = dfstrat[dfstrat$country == unique(dfstrat$country)[1] & dfstrat$temps == years,]$valeur
  temp_df = data.frame(valeur=newrow, country=unique(dfstrat$country)[1], temps = years, index=unique(df.melt[df.melt$country==unique(dfstrat$country)[1],]$index), 
weight = dfstrat[dfstrat$country == unique(dfstrat$country)[1] & dfstrat$temps == years, ]$weight)
  dfstratfin = rbind(dfstratfin, temp_df)
  
  for (countries in unique(dfstrat$country)[2:nb_country_plot]){ #on garde ceux les plus élevés
    newrow = newrow + dfstrat[dfstrat$country == countries & dfstrat$temps == years, ]$valeur
    
    temp_df = data.frame(valeur=newrow, index=unique(df.melt[df.melt$country==countries,]$index), country=countries, 
weight = dfstrat[dfstrat$country == countries & dfstrat$temps == years,]$weight, temps = years)
    dfstratfin = rbind(dfstratfin, temp_df)
  }}



#On rajoute l'unité synthétique
for (years in date0:date3){
  newrow = dfstrat[dfstrat$country == unique(dfstrat$country)[1] & dfstrat$temps == years,]$valeur
  for (countries in unique(dfstrat$country)[2:length(unique(dfstrat$country))]){ #on rajoute les petites valeurs
    newrow = newrow + dfstrat[dfstrat$country == countries & dfstrat$temps == years, ]$valeur}
  temp_df = data.frame(valeur=newrow, index=0, country='Synthetic Unit', weight = 1, temps = years)
  dfstratfin = rbind(dfstratfin, temp_df)
}
dfstratfin = dfstratfin[order(dfstratfin$weight, decreasing=TRUE),]


#On ajoute la vraie valeur pour la France
for (years in date0:date3){
  temp_df = data.frame(valeur=df.melt[df.melt$country == pays & df.melt$temps == years,][,valeur], country= pays, temps=years, weight = 1, index = ind_pays)

  dfstratfin=rbind(dfstratfin, temp_df)}

colorplot = ggplot(data=dfstratfin,aes(x=temps, y=valeur))+
  geom_line(aes(group=country, col=country))

print(colorplot)
}

```
```{r}
df_stratTitouan()
```





###Fonction synth sans avoir à tout mettre
```{r}
simpleSynth=function(liste_value, df.melt=df.melted_i, pays='FRA', dateplot=1997:2016, dateopti=2000:2012){


special=list()
###création d'une liste pour les prédicteurs spéciaux
year=seq(dateopti[1],dateopti[length(dateopti)],2)


ind_pays=which(df.melt$country==pays)[1]
#Version avec les différentes dates version Synth
dataprep.out=Synth::dataprep(foo=df.melt,
                predictors.op = "mean", #La manière dont on les prend en compte
                time.predictors.prior = 2012, #période pré-traitement
                special.predictors = gen_liste(liste_value),
                dependent=liste_value[1], #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = ind_pays, # n° de l'unité traitée
                controls.identifier = c(1:max(df.melt$index))[-ind_pays], #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = dateopti, #période de maximisation
                time.plot=dateplot) #période tracée sur le graphique

synth.out <- Synth::synth(data.prep.obj = dataprep.out, method = "BFGS")

Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = liste_value[1], Xlab = "year", Legend = c(pays, paste('synthetic',pays)), Legend.position = "bottomright")
synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
poids1=synth.tables$tab.w%>%rename(poids=w.weights)
print(poids1)
df_stratTitouan(df.melt=df.melt, liste_value=liste_value, pays=pays, date0=dateplot[1], date1=dateopti[1], date2=dateopti[length(dateopti)], date3=dateplot[length(dateplot)], synth.out= synth.out, dataprep.out=dataprep.out, synthdon=TRUE)}
```
```{r}
simpleSynth(c('evol_VA', 'W', 'Productivite','Invest_pop'), df.melt=df.evol.filtered, dateplot=2000:2016)
```



###Une fonction pour afficher le contrôle synthétique avec la fonction de base, celle en s'apesantissant sur les poids et celle en calculant les poids à part
```{r}
tripleSynth=function(liste_value, df.melt, pays='FRA', date0=1999, date1=2000, date2=2012, date3=2016){

#calcul des poids :
opti=calculpoids(df.melt, liste_value[1], pays)
special=list()
###création d'une liste pour les prédicteurs spéciaux
year=seq(date0,date2,2)

ind_pays=as.numeric(df.melt[df.melt$country==pays,]$index[1])
#Version avec les différentes dates version Synth
dataprep.out=Synth::dataprep(foo=df.melt,
                predictors.op = "mean", #La manière dont on les prend en compte
                time.predictors.prior = 2012, #période pré-traitement
                special.predictors = gen_liste(liste_value),
                dependent=liste_value[1], #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = ind_pays, # n° de l'unité traitée
                controls.identifier = c(1:max(df.melt$index))[-ind_pays], #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = date1:date2, #période de maximisation
                time.plot=date0:date3) #période tracée sur le graphique

synth.out <- Synth::synth(data.prep.obj = dataprep.out, method = "BFGS")

Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = liste_value[1], Xlab = "year", Legend = c(pays, paste('synthetic',pays)), Legend.position = "bottomright", Main='avec toutes les dates')
synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
poids1=synth.tables$tab.w%>%rename(poids=w.weights)
#Version sans les différentes dates
dataprep.out=Synth::dataprep(foo=df.melt,
                predictors = liste_value, #Les variables utilisées en contrôle
                predictors.op = "mean", #La manière dont on les prend en compte
                time.predictors.prior = 2012, #période pré-traitement
                dependent=liste_value[1], #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = ind_pays, # n° de l'unité traitée
                controls.identifier = c(1:max(df.melt$index))[-ind_pays], #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = date1:date2, #période de maximisation
                time.plot=date0:date3) #période tracée sur le graphique
synth.out <- Synth::synth(data.prep.obj = dataprep.out, method = "BFGS")

Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = liste_value[1], Xlab = "year", Legend = c(pays, paste('synthetic',pays)), Legend.position = "bottomright", Main='avec la moyenne')
 
synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,synth.res = synth.out)
pds=synth.tables$tab.w%>%inner_join(poids1, on='unit.names')
pds$poids_main=opti
print(pds)

#version a la mano
synth.out$solution.w=opti#on remplace les poids par ceux calculés

Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = liste_value[1], Xlab = "year", Legend = c(pays, paste('synthetic',pays)), Legend.position = "bottomright", Main ='à la main')
}
tripleSynth(liste_value=c("VA_ind"), df.melted_i, pays='GBR', date1=2005)
```




```{r}
tripleSynth(liste_value=c("Emp_ind","VA_ind", "Productivite", "export", "W", "taux_imposition"), df.fisc, pays='FRA', date0=2000, date1=2005)
```
```{r}

tripleSynth(liste_value=c("evol_Emp","evol_VA"), df.evol, pays='FRA', date0=2000, date1=2005)
```














```{r}
variable=c("Invest_pop","W","Productivite","export","Part_lab","chomage","Emp_ind")
df.melted_i$Emploi_scale<-scale(df.melted_i$Emploi)

indice_sans=function(base)
{
 df_exo=filter(df.melted_i, ! country %in% c('DEU','ITA','PRT'))
 a=unique(df_exo$index) 
 a=a[a!=11]
 print(a)
 return(a)  
}
d=c(1:10,12:24)[-indice_sans(df.melted_i)]
b=c(1:10,12:24)

dataprep.out=Synth::dataprep(foo=df.melted_i,
                predictors = (variable),#Les variables utilisées en contrôle
                time.predictors.prior = c(2000:2011), #période pré-traitement
                special.predictors = list(list("VA_ind",seq(2000,2007),"mean"),
                                          list("VA_ind",seq(2010,2011),"mean")), 
                dependent="VA_ind", #variable d'intérêt
                unit.variable = "index", #colonne avec les index des pays
                unit.names.variable="country", #nom des pays
                treatment.identifier = 11, # n° de l'unité traitée
                controls.identifier = b[-d], #n° des unités de contrôles
                time.variable = 'temps', #colonne où se trouve la date
                time.optimize.ssr = 2000:2010, #période de maximisation
                time.plot=1999:2016,
                ) #période tracée sur le graphique


synth.out <- synthBound(data.prep.obj = dataprep.out, quadopt = 'ipop')

Synth::path.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = "Emploi indice 100", Xlab = "year", Legend = c("France", "synthetic France"), Legend.position = "bottomleft")
Synth::gaps.plot(synth.res = synth.out, dataprep.res = dataprep.out,
Ylab = "Emploi indice 100", Xlab = "year")

synth.tables <- Synth::synth.tab(dataprep.res = dataprep.out,
 synth.res = synth.out)

synth.tables$tab.w%>%arrange(-synth.tables$tab.w$w.weight)
synth.tables$tab.v
```



```{r}
write.table(df.melted_i, "donnee.csv",sep=",")
read.csv("donnee.csv")
```


```{r}
Synth::gaps.plot(synth.res = synth.out, dataprep.res = dataprep.out,Ylab = "PIB", Xlab = "year")

valeur=dataprep.out$Y1-dataprep.out$Y0%*%synth.out$solution.w
```


```{r}
numinit=2
numend=20
d=c(1:24)
variable=c("Invest_pop","W","VA_ind","Productivite","export","Part_lab","chomage")
df.melted_i
for (i in c(numinit:numend))
{
  #On calcul le contrôle synthétique avec l'unité i traitée
  print(i)
  ind_ctrl=c(1:(i-1),(i+1):21)
  dataprep_i.out=Synth::dataprep(foo=df.melted_i,
                predictors = (variable),
                special.predictors = list(list("VA_ind",seq(2000,2005),"mean"),
                                          list("VA_ind",seq(2007,2008),"mean"),
                                          list("VA_ind",seq(2010,2011),"mean")),
                time.predictors.prior = 2000:2011, 
                dependent="VA_ind", 
                unit.variable = "index", 
                unit.names.variable="country", 
                treatment.identifier = i, 
                controls.identifier =   ind_ctrl, 
                time.variable = 'temps', 
                time.optimize.ssr = 2000:2010, #
                time.plot=1999:2016,) 

quiet(synth.out_i <- Synth::synth(data.prep.obj = dataprep_i.out, method = "BFGS")) #On calcul les résultats sans les afficher

#Synth::gaps.plot(synth.res = synth.out_i, dataprep.res = dataprep_i.out,Ylab = "PIB", Xlab = "year")

  if (i==numinit) { #On initialise une base avec en ligne le noms de l'unité puis l'effet du traitement pour chaque année (le numéro de l'année n'est pas dans la base)
    noms=c(rep(as.character(i),nrow(dataprep_i.out$Y0)))
    valeur=dataprep_i.out$Y1-dataprep_i.out$Y0%*%synth.out_i$solution.w 
    df=data.frame(noms,valeur)
    couleur=c("red") #ca c'est pour gérer les couleurs 
  }
    
  if (i>numinit) {  #On complète la base
    noms=c(noms,rep(as.character(i),nrow(dataprep_i.out$Y0)))
    valeur=c(valeur,dataprep_i.out$Y1-dataprep_i.out$Y0%*%synth.out_i$solution.w)
    df=data.frame(noms,valeur)
    if (i==10){couleur=c(couleur,"blue")}
    else{couleur=c(couleur,"red")}}

}
    
df
df2=data.frame(x=(rep((1999:2016),length(unique(noms)))),df)
df_effe=filter(df2,x>2008)
```


```{r}
df2$noms=as.numeric(df2$nom)
df_RMPSE=df2%>%rename(index=noms)%>%left_join(df.melted_i[c('country','index')][1:24,])
df_RMPSE$valeur_sq=df_RMPSE$valeur*df_RMPSE$valeur
RMPSE_res=group_by(df_RMPSE,index)%>% summarise(RMPSE_tot = mean(valeur_sq),country=unique(country))
RMPSE_res$RMPSE_trait=(group_by(filter(df_RMPSE,x>2012),index)%>% summarise(mean = mean(valeur_sq),country=unique(country)))$mean
RMPSE_res$RMPSE_prev=(group_by(filter(df_RMPSE,x<2013),index)%>% summarise(mean = mean(valeur_sq),country=unique(country)))$mean
```


```{r}
RMPSE_res
RMPSE_res%>%mutate(RMPSE_tot=sqrt(RMPSE_res$RMPSE_tot),
                   RMPSE_trait=sqrt(RMPSE_res$RMPSE_trait),
                   RMPSE_res=sqrt(RMPSE_res$RMPSE_prev)) %>%select(country,RMPSE_tot,RMPSE_trait,RMPSE_prev)
RMPSE_res$RMPSE_ajuste=RMPSE_res$RMPSE_trait/RMPSE_res$RMPSE_prev
RMPSE_res%>%arrange(-RMPSE_res$RMPSE_ajuste)
```


```{r}
df2$noms=as.character(df2$noms)
qplot(x=x, y=valeur,color=noms, data=df2, geom="line",ylim=c(-20,25))+ xlab("Date") +scale_color_manual(breaks = unique(df2$noms),                        values=couleur)
df.melted_i
```

###Tentative en normalisant###
```{r}
SynthNorm=function(value, df.melt=df.melted,pays='FRA', date1=1999, date2=2013, date0=1997, date3=2016){
pays_autres=length(c(table(df.melt$country)))-1
Xt=subset(df.melt, country==pays& temps<date2 &temps>date1)
Xt=Xt[,value]
Xtnorm=scale(Xt)[1:(date2-date1-1)]#on récupère les données pour la France

Xn=subset(df.melt, temps<date2 &temps>date1)
Xn=Xn[Xn$country!=pays,]#on récupère les autres
Xn=matrix(Xn[,value], nrow=pays_autres)
Xnorm=t(scale(t(Xn)))#on normalise les données par pays(la normalisation se fait par colonne)

#calcul des poids
x0=c(rep(1/pays_autres,pays_autres))#on crée un vecteur de départ au hasard
f=function(x){norm(Xtnorm-x%*%Xnorm, type="2")}
n=function(x){x}#contrainte de positivité des poids (hin(x)>0)
q=function(x){sum(x)-1}#fonction de la contrainte d'égalité (heq(x)=0)
opti=quiet(auglag(x0, f, hin=n, heq=q))#fonction d'optimisation sur toute la période
opti$par=matrix(opti$par, ncol=1)
#plot
Xt=subset(df.melt, country==pays)
Xt=Xt[,value]
Xn=subset(df.melt)
Xn=Xn[Xn$country!=pays,]#on récupère les autres
Xn=matrix(Xn[,value], nrow=pays_autres)
Xnorm=t(scale(t(Xn)))#on normalise les données par pays(la normalisation se fait par colonne)
FrSynth=t(opti$par)%*%Xnorm*sd(c(Xt))+mean(c(Xt))
df=data.frame(c(FrSynth, Xt))
df$X=c(rep(2,nrow(df)/2),rep(3,nrow(df)/2))
df$dat=c(date0:date3, date0:date3)
ggplot(data=df,aes(x=dat, c.FrSynth..Xt., col=X)) +
  geom_line(aes(group=X))
}
SynthNorm('VA')
```

```{r}
SynthNorm('chomage',df.melt=df.melted_i)
```

```{r}
date1=1999
date2=2013
date0=1997
date3=2016
value='VA'
df.melt=df.melted
pays_autres=length(c(table(df.melt$country)))-1
Xt=subset(df.melt, country==pays& temps<date2 &temps>date1)
Xt=Xt[,value]
Xtnorm=scale(Xt)[1:(date2-date1-1)]#on récupère les données pour la France

Xn=subset(df.melt, temps<date2 &temps>date1)
Xn=Xn[Xn$country!=pays,]#on récupère les autres
Xn=matrix(Xn[,value], nrow=pays_autres)
Xnorm=t(scale(t(Xn)))#on normalise les données par pays(la normalisation se fait par colonne)
```

#calcul des poids
```{r}
x0=c(rep(1/pays_autres,pays_autres))#on crée un vecteur de départ au hasard
f=function(x){norm(Xtnorm-x%*%Xnorm, type="2")}
n=function(x){x}#contrainte de positivité des poids (hin(x)>0)
q=function(x){sum(x)-1}#fonction de la contrainte d'égalité (heq(x)=0)
opti=quiet(auglag(x0, f, hin=n, heq=q))#fonction d'optimisation sur toute la période
opti$par=matrix(opti$par, ncol=1)
```

#plot des résultats
```{r}
Xt=subset(df.melted, country==pays)
Xt=Xt[,value]
Xn=subset(df.melted)
Xn=Xn[Xn$country!=pays,]#on récupère les autres
Xn=matrix(Xn[,value], nrow=pays_autres)
Xnorm=t(scale(t(Xn)))#on normalise les données par pays(la normalisation se fait par colonne)
FrSynth=t(opti$par)%*%Xnorm*sd(c(Xt))+mean(c(Xt))
plot(c(FrSynth))
plot(Xt)
df=data.frame(c(FrSynth, Xt))
df$X=c(rep(2,nrow(df)/2),rep(3,nrow(df)/2))
df$dat=c(1997:2016, 1997:2016)
ggplot(data=df,aes(x=dat, c.FrSynth..Xt., col=X)) +
  geom_line(aes(group=X))

```

```{r}
synthBound <-
function(           data.prep.obj = NULL,
                      X1 = NULL,
                      X0 = NULL,
                      Z0 = NULL,
                      Z1 = NULL,
                      custom.v = NULL,
                      optimxmethod = c("Nelder-Mead","BFGS"),
                      genoud = FALSE,
                      quadopt = "ipop",
                      Margin.ipop = 0.0005,
                      Sigf.ipop = 5,
                      Bound.ipop = 10,
                      verbose = FALSE,
                       ...
                      )
  { 
    # Retrieve dataprep objects
    if(is.null(data.prep.obj) == FALSE)
      {
        cat("\nX1, X0, Z1, Z0 all come directly from dataprep object.\n\n")
        X1 <- data.prep.obj$X1
        Z1 <- data.prep.obj$Z1
        X0 <- data.prep.obj$X0
        Z0 <- data.prep.obj$Z0
      } else {
        cat("X1,X0,Z1,Z0 were individually input (not dataprep object.)\n\n")
      }

     # routine checks
     store <- list(X1=X1,X0=X0,Z1=Z1,Z0=Z0)
     for(i in 1:4){
     if(is.null(store[[i]]))
      {stop(paste("\n",names(store)[i],"is missing \n"))}
     if(sum(is.na(store[[i]]))>0)
      {stop(paste("\n NAs in",names(store)[i],"\n"))}
     if(is.matrix(store[[i]]) == FALSE)
      {stop(paste("\n",names(store)[i],"is not a matrix object\n"))}
     }

    # geometry checks
    if(ncol(X1)!=1){stop("\n Please specify only one treated unit: X1 has to have ncol= 1")}
    if(ncol(Z1)!=1){stop("\n Please specify only one treated unit: Z1 has to have ncol= 1")}

    if(ncol(X0)<2){stop("\n Please specify at least two control units: X0 has to have ncol >= 2 ")}
    if(ncol(Z0)<2){stop("\n Please specify only one treated unit: Z0 has to have ncol >= 2")}

    if(nrow(Z0)!=nrow(Z1)){stop("\n Different number of periods for treated and controls: nrow(Z0) unequal nrow(Z1)")}
    if(nrow(X0)!=nrow(X1)){stop("\n Different number of predictors for treated and controls: nrow(X0) unequal nrow(X1)")}

    if(nrow(X0)==0){stop("No predictors specified. Please specify at least on predictor")}
    if(nrow(Z0)==0){stop("No periods specified for Z1 and Z0. Please specify at least on period")}

    if(0 %in% apply(X0,1,sd))
     {stop("\n At least one predictor in X0 has no variation across control units. Please remove this predictor.")}

    # collinearity check
   # check <- try(solve(t(X0)%*%X0),silent=TRUE)
   # if(class(check)=="try-error")
   #  {stop("\n Some of the predictors in X0 are collinear (t(X0)%*%X0) not invertible")}

    # Normalize X
    nvarsV <- dim(X0)[1]
    big.dataframe <- cbind(X0, X1)
    divisor <- sqrt(apply(big.dataframe, 1, var))
    scaled.matrix <-
      t(t(big.dataframe) %*% ( 1/(divisor) *
                              diag(rep(dim(big.dataframe)[1], 1)) ))

    X0.scaled <- scaled.matrix[,c(1:(dim(X0)[2]))]
    if(is.vector(X0.scaled)==TRUE)
     {X0.scaled <- t(as.matrix(X0.scaled))}
    X1.scaled <- scaled.matrix[,dim(scaled.matrix)[2]]


    # check if custom v weights are supplied or
    # if only on predictor is specified,
    # we jump to quadratic optimization over W weights
    # if not start optimization over V and W
    if(is.null(custom.v) & nrow(X0) != 1)
      {
      
      # two attemps for best V are made: 
      # equal weights and regression based starting values 
     cat("\n****************",
         "\n searching for synthetic control unit  \n","\n"
        )
      
      if(genoud == TRUE) # if user wants genoud as well
      {
      require(rgenoud)                # we run genoud first
      cat("\n****************",
         "\n genoud() requested for optimization\n","\n"
        )

      rgV.genoud <- genoud(
                             fn.V, 
                             nvarsV, 
                             X0.scaled = X0.scaled,
                             X1.scaled = X1.scaled,
                             Z0 = Z0,
                             Z1 = Z1,
                             quadopt = quadopt,
                             margin.ipop = Margin.ipop,
                             sigf.ipop = Sigf.ipop,
                             bound.ipop = Bound.ipop
                             )
      SV1 <- rgV.genoud$par  # and use these as starting values

      cat("\n****************",
         "\n genoud() finished, now running local optimization using optim()\n","\n"
        )

      } else {
      # if we don't use genoud first: set of starting values: equal weights
      SV1 <- rep(1/nvarsV,nvarsV)
      }
      
      # now we run optimization
      all.methods <- FALSE
      if(sum(optimxmethod %in% c("All"))==1){ all.methods <- TRUE }
     rgV.optim.1 <- optimx(par=SV1, fn=fn.V,
                             gr=NULL, hess=NULL, 
                             method=optimxmethod, itnmax=NULL, hessian=FALSE,
                             control=list(kkt=FALSE,
                                          starttests=FALSE,
                                          dowarn=FALSE,
                                          all.methods=all.methods),
                             X0.scaled = X0.scaled,
                             X1.scaled = X1.scaled,
                             Z0 = Z0,
                             Z1 = Z1,
                             quadopt = quadopt,
                             margin.ipop = Margin.ipop,
                             sigf.ipop = Sigf.ipop,
                             bound.ipop = Bound.ipop
                            )
      # get minimum
      if(verbose==TRUE){print(rgV.optim.1)}
      rgV.optim.1 <- collect.optimx(rgV.optim.1,"min")
      
      # second set of starting values: regression method 
      # will sometimes not work because of collinear Xs
      # so it's wrapped in a try command
      Xall <- cbind(X1.scaled,X0.scaled)
      Xall <- cbind(rep(1,ncol(Xall)),t(Xall))
      Zall <- cbind(Z1,Z0)
      Beta <- try(solve(t(Xall)%*%Xall)%*%t(Xall)%*%t(Zall),silent=TRUE)
      
      # if inverses did not work, we
      # stick with first results    
      if(inherits(Beta,"try-error")) 
       {
        rgV.optim <- rgV.optim.1
       } else {
      # otherwise we run a second optimization with regression starting values
        Beta <- Beta[-1,]
        V    <- Beta%*%t(Beta)
        SV2  <- diag(V)
        SV2 <- SV2 / sum(SV2)
  
      rgV.optim.2 <- optimx(par=SV2, fn=fn.V,
                             gr=NULL, hess=NULL, 
                             method=optimxmethod, itnmax=NULL, hessian=FALSE,
                             control=list(kkt=FALSE,
                                          starttests=FALSE,
                                          dowarn=FALSE,
                                          all.methods=all.methods),
                             X0.scaled = X0.scaled,
                             X1.scaled = X1.scaled,
                             Z0 = Z0,
                             Z1 = Z1,
                             quadopt = quadopt,
                             margin.ipop = Margin.ipop,
                             sigf.ipop = Sigf.ipop,
                             bound.ipop = Bound.ipop
                            )
      # get minimum
      if(verbose==TRUE){print(rgV.optim.2)}
      rgV.optim.2 <- collect.optimx(rgV.optim.2,"min")
  
      # ouput
      if(verbose == TRUE){
      cat("\n Equal weight loss is:",rgV.optim.1$value,"\n")
      cat("\n Regression Loss is:",rgV.optim.2$value,"\n")
      }       
      # and keep the better optim results    
      if(rgV.optim.1$value < rgV.optim.2$value) 
       {
        rgV.optim <- rgV.optim.1
       } else {
        rgV.optim <- rgV.optim.2
       }
      } # close if statement for second regression based optimization attempt
     
      # final V weights from optimization
      solution.v   <- abs(rgV.optim$par)/sum(abs(rgV.optim$par))
     } else { # jump here if only optimize over W


     cat("\n****************",
         "\n optimization over w weights: computing synthtic control unit \n","\n\n"
         )

     if(nrow(X0)==1)
      {
       custom.v <- 1 # this is the case where only one predictor is specified V is the identity matrix
      } else {
     # if user choose to supply v manually:
     cat("\n****************",
         "\n v weights supplied manually: computing synthtic control unit \n","\n\n"
         )
     if(length(custom.v) != nvarsV){stop("custom.V misspecified: length(custom.V) != nrow(X1)")}
     if(mode(custom.v)!="numeric"){stop("custom.V must be numeric")}
      }

    # enter solution.V
    rgV.optim  <- NULL
    solution.v <- abs(custom.v)/sum(custom.v)
    
  } # close else statment for by-passing V optimization


    # last step: now recover solution.w
    V <- diag(x=as.numeric(solution.v),nrow=nvarsV,ncol=nvarsV)
    H <- t(X0.scaled) %*% V %*% (X0.scaled)
    a <- X1.scaled
    c <- -1*c(t(a) %*% V %*% (X0.scaled) )
    A <- t(rep(1, length(c)))
    b <- 1
    l <- rep(0, length(c))
    u <- rep(0.25, length(c))
    r <- 0

     if(quadopt=="ipop"){
    res <- ipop(c = c, H = H, A = A, b = b, l = l, u = u, r = r,
                margin = Margin.ipop, maxiter = 1000, sigf = Sigf.ipop, bound = Bound.ipop)
    solution.w <- as.matrix(primal(res))
    } else {
    # LowRankQP
     if(quadopt=="LowRankQP"){
      res <- LowRankQP(Vmat=H,dvec=c,Amat=A,bvec=1,uvec=rep(1,length(c)),method="LU")
      solution.w <- as.matrix(res$alpha)
     } 
    }

    rownames(solution.w) <- colnames(X0)
    colnames(solution.w) <- "w.weight"
    names(solution.v) <- rownames(X0)

    loss.w <- t(X1.scaled - X0.scaled %*% solution.w) %*%
      V %*% (X1.scaled - X0.scaled %*% solution.w)

    loss.v <-
      t(Z1 - Z0 %*% as.matrix(solution.w)) %*%
        (Z1 - Z0 %*% as.matrix(solution.w)) 
    loss.v <- loss.v/nrow(Z0)      
 
    # produce viewable output
    cat("\n****************",
        "\n****************",
        "\n****************",
        "\n\nMSPE (LOSS V):", loss.v,
  #      "\n\nLOSS (W):", loss.w,
        "\n\nsolution.v:\n", round(as.numeric(solution.v), 10),
        "\n\nsolution.w:\n", round(as.numeric(solution.w), 10),
        "\n\n"
        )
        
    optimize.out <- list(
                         solution.v = solution.v,
                         solution.w = solution.w,
                         loss.v = loss.v,
                         loss.w = loss.w,
                         custom.v = custom.v,
                         rgV.optim = rgV.optim
                         )

    return(invisible(optimize.out))

  }
```